<commit_before>/*
 * 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 <vector>

#include "CubemapIBL.h"

#include "math/mat3.h"
#include <math/scalar.h>
#include <utils/JobSystem.h>

#include "Cubemap.h"
#include "CubemapUtils.h"
#include "ProgressUpdater.h"
#include "utilities.h"

using namespace math;
using namespace image;
using namespace utils;

extern bool g_quiet;

static double pow5(double x) {
    return (x*x)*(x*x)*x;
}

static double3 hemisphereImportanceSampleDggx(double2 u, double a) {
    const double phi = 2 * M_PI * u.x;
    // NOTE: (aa-1) == (a-1)(a+1) produces better fp accuracy
    const double cosTheta2 = (1 - u.y) / (1 + (a + 1) * ((a - 1) * u.y));
    const double cosTheta = std::sqrt(cosTheta2);
    const double sinTheta = std::sqrt(1 - cosTheta2);
    return { sinTheta * std::cos(phi), sinTheta * std::sin(phi), cosTheta };
}

static double3 __UNUSED hemisphereCosSample(double2 u) {
    const double phi = 2 * M_PI * u.x;
    const double cosTheta2 = 1 - u.y;
    const double cosTheta = std::sqrt(cosTheta2);
    const double sinTheta = std::sqrt(1 - cosTheta2);
    return { sinTheta * std::cos(phi), sinTheta * std::sin(phi), cosTheta };
}

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

static double __UNUSED DistributionAshikhmin(double NoH, double linearRoughness) {
    double a = linearRoughness;
    double a2 = a * a;
    double cos2h = NoH * NoH;
    double sin2h = 1 - cos2h;
    double sin4h = sin2h * sin2h;
    return 1 / (M_PI * (1 + 4 * a2)) * (sin4h + 4 * std::exp(-cos2h / (a2 * sin2h)));
}

static double __UNUSED DistributionCharlie(double NoH, double linearRoughness) {
    // Estevez and Kulla 2017, "Production Friendly Microfacet Sheen BRDF"
    double a = linearRoughness;
    double invAlpha = 1 / a;
    double cos2h = NoH * NoH;
    double sin2h = 1 - cos2h;
    return (2 + invAlpha) * std::pow(sin2h, invAlpha * 0.5) / (2 * M_PI);
}

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

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

static double __UNUSED VisibilityAshikhmin(double NoV, double NoL, double a) {
    return 1 / (4 * (NoL + NoV - NoL * NoV));
}

void CubemapIBL::roughnessFilter(Cubemap& dst,
        const std::vector<Cubemap>& levels, double linearRoughness, size_t maxNumSamples)
{
    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 = float((4 * M_PI) / (6 * dim0 * dim0));

    // be careful w/ the size of this structure, the smaller the better
    struct CacheEntry {
        double3 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 #
    double 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, sample = 0 ; sampleIndex < maxNumSamples; sampleIndex++) {
        /*
         *       (sampling)
         *            L         H (never calculated below)
         *            .\       /.
         *            . \     / .
         *            .  \   /  .
         *            .   \ /   .
         *         ---|----o----|-------> n
         *    cos(2*theta)    cos(theta)
         *       = n.L           = n.H
         *
         * Note: NoH == LoH
         * (H is the half-angle between L and V, and V == N)
         *
         */

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

        // Importance sampling GGX - Trowbridge-Reitz
        // This corresponds to integrating Dggx()*cos(theta) over the hemisphere
        const double3 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 double3 N = {0, 0, 1};
        const double3 V = N;
        const double3 L = 2 * dot(H, V) * H - V;
        const double NoL = dot(N, L);
        const double NoH = dot(N, H);
        const double NoH2 = NoH * NoH;
        const double NoV = dot(N, V);
        const double LoH = dot(L, H);
#else
        const double NoV = 1;
        const double NoH = H.z;
        const double NoH2 = H.z*H.z;
        const double NoL = 2*NoH2 - 1;
        const double3 L(2*NoH*H.x, 2*NoH*H.y, NoL);
        const double LoH = dot(L, H);
#endif

        // pre-filtered importance sampling
        // see: "Real-time Shading with Filtered Importance Sampling", <NAME>
        // see: "GPU-Based Importance Sampling, GPU Gems 3", <NAME>

        // K is a LOD bias that allows a bit of overlapping between samples
        // log4(K)=1 empirically works well with box-filtered mipmaps
        constexpr float K = 4;

        // OmegaS is is the solid-angle of an important sample (i.e. how much surface
        // of the sphere it represents). It obviously is function of the PDF.
        const double pdf = DistributionGGX(NoH, linearRoughness) / 4;
        const double omegaS = 1 / (numSamples * pdf);

        // The LOD is given by: max[ log4(Os/Op) + K, 0 ]
        const double l = float(log4(omegaS) - log4(omegaP) + log4(K));
        const float mipLevel = clamp(float(l), 0.0f, maxLevelf);

        if (NoL > 0) {
            const double V = Visibility(NoV, NoL, linearRoughness);
            const double F = Fresnel(1, 0, LoH);
            const float brdf_NoL = float(F * V * 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 - l0;
            cache.push_back({ L, brdf_NoL, lerp, l0, l1 });
            sample++;
        }
    }

    std::for_each(cache.begin(), cache.end(), [weight](CacheEntry& entry){
        entry.brdf_NoL /= 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;
    });

    ProgressUpdater updater(1);
    if (!g_quiet) {
        updater.start();
    }<commit_msg>Fold computations inside the branch<commit_after>/*
 * 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 <vector>

#include "CubemapIBL.h"

#include "math/mat3.h"
#include <math/scalar.h>
#include <utils/JobSystem.h>

#include "Cubemap.h"
#include "CubemapUtils.h"
#include "ProgressUpdater.h"
#include "utilities.h"

using namespace math;
using namespace image;
using namespace utils;

extern bool g_quiet;

static double pow5(double x) {
    return (x*x)*(x*x)*x;
}

static double3 hemisphereImportanceSampleDggx(double2 u, double a) {
    const double phi = 2 * M_PI * u.x;
    // NOTE: (aa-1) == (a-1)(a+1) produces better fp accuracy
    const double cosTheta2 = (1 - u.y) / (1 + (a + 1) * ((a - 1) * u.y));
    const double cosTheta = std::sqrt(cosTheta2);
    const double sinTheta = std::sqrt(1 - cosTheta2);
    return { sinTheta * std::cos(phi), sinTheta * std::sin(phi), cosTheta };
}

static double3 __UNUSED hemisphereCosSample(double2 u) {
    const double phi = 2 * M_PI * u.x;
    const double cosTheta2 = 1 - u.y;
    const double cosTheta = std::sqrt(cosTheta2);
    const double sinTheta = std::sqrt(1 - cosTheta2);
    return { sinTheta * std::cos(phi), sinTheta * std::sin(phi), cosTheta };
}

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

static double __UNUSED DistributionAshikhmin(double NoH, double linearRoughness) {
    double a = linearRoughness;
    double a2 = a * a;
    double cos2h = NoH * NoH;
    double sin2h = 1 - cos2h;
    double sin4h = sin2h * sin2h;
    return 1 / (M_PI * (1 + 4 * a2)) * (sin4h + 4 * std::exp(-cos2h / (a2 * sin2h)));
}

static double __UNUSED DistributionCharlie(double NoH, double linearRoughness) {
    // Estevez and Kulla 2017, "Production Friendly Microfacet Sheen BRDF"
    double a = linearRoughness;
    double invAlpha = 1 / a;
    double cos2h = NoH * NoH;
    double sin2h = 1 - cos2h;
    return (2 + invAlpha) * std::pow(sin2h, invAlpha * 0.5) / (2 * M_PI);
}

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

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

static double __UNUSED VisibilityAshikhmin(double NoV, double NoL, double a) {
    return 1 / (4 * (NoL + NoV - NoL * NoV));
}

void CubemapIBL::roughnessFilter(Cubemap& dst,
        const std::vector<Cubemap>& levels, double linearRoughness, size_t maxNumSamples)
{
    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 = float((4 * M_PI) / (6 * dim0 * dim0));

    // be careful w/ the size of this structure, the smaller the better
    struct CacheEntry {
        double3 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 #
    double 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, sample = 0 ; sampleIndex < maxNumSamples; sampleIndex++) {
        /*
         *       (sampling)
         *            L         H (never calculated below)
         *            .\       /.
         *            . \     / .
         *            .  \   /  .
         *            .   \ /   .
         *         ---|----o----|-------> n
         *    cos(2*theta)    cos(theta)
         *       = n.L           = n.H
         *
         * Note: NoH == LoH
         * (H is the half-angle between L and V, and V == N)
         *
         */

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

        // Importance sampling GGX - Trowbridge-Reitz
        // This corresponds to integrating Dggx()*cos(theta) over the hemisphere
        const double3 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 double3 N = {0, 0, 1};
        const double3 V = N;
        const double3 L = 2 * dot(H, V) * H - V;
        const double NoL = dot(N, L);
        const double NoH = dot(N, H);
        const double NoH2 = NoH * NoH;
        const double NoV = dot(N, V);
        const double LoH = dot(L, H);
#else
        const double NoV = 1;
        const double NoH = H.z;
        const double NoH2 = H.z*H.z;
        const double NoL = 2*NoH2 - 1;
        const double3 L(2*NoH*H.x, 2*NoH*H.y, NoL);
        const double LoH = dot(L, H);
#endif

        if (NoL > 0) {
            // pre-filtered importance sampling
            // see: "Real-time Shading with Filtered Importance Sampling", <NAME>
            // see: "GPU-Based Importance Sampling, GPU Gems 3", <NAME>

            // K is a LOD bias that allows a bit of overlapping between samples
            // log4(K)=1 empirically works well with box-filtered mipmaps
            constexpr float K = 4;

            // OmegaS is is the solid-angle of an important sample (i.e. how much surface
            // of the sphere it represents). It obviously is function of the PDF.
            const double pdf = DistributionGGX(NoH, linearRoughness) / 4;
            const double omegaS = 1 / (numSamples * pdf);

            // The LOD is given by: max[ log4(Os/Op) + K, 0 ]
            const double l = float(log4(omegaS) - log4(omegaP) + log4(K));
            const float mipLevel = clamp(float(l), 0.0f, maxLevelf);

            const double V = Visibility(NoV, NoL, linearRoughness);
            const double F = Fresnel(1, 0, LoH);
            const float brdf_NoL = float(F * V * 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 - l0;

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

            sample++;
        }
    }

    std::for_each(cache.begin(), cache.end(), [weight](CacheEntry& entry){
        entry.brdf_NoL /= 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;
    });

    ProgressUpdater updater(1);
    if (!g_quiet) {
        updater.start();
    }
{'fraction_non_alphanumeric': 0.11803696622052262, 'fraction_numerical': 0.020395156150414276, 'mean_word_length': 2.299894847528917, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3843377', 'n_tokens_mistral': 5389, 'n_tokens_neox': 5126, 'n_words': 1690}