diff options
Diffstat (limited to 'intern/cycles/app/cycles_precompute.cpp')
| -rw-r--r-- | intern/cycles/app/cycles_precompute.cpp | 400 |
1 files changed, 400 insertions, 0 deletions
diff --git a/intern/cycles/app/cycles_precompute.cpp b/intern/cycles/app/cycles_precompute.cpp new file mode 100644 index 00000000000..6f9d9601b1e --- /dev/null +++ b/intern/cycles/app/cycles_precompute.cpp @@ -0,0 +1,400 @@ +#include "util/math.h" +#include "util/string.h" +#include "util/system.h" + +#include "util/hash.h" +#include "util/task.h" + +#include "kernel/device/cpu/compat.h" +#include "kernel/device/cpu/globals.h" + +#include "kernel/sample/lcg.h" +#include "kernel/sample/mapping.h" + +#include "kernel/closure/bsdf_microfacet.h" +#include "kernel/closure/bsdf_microfacet_glass.h" +#include "kernel/closure/bsdf_principled_sheen.h" + +#include <iostream> + +CCL_NAMESPACE_BEGIN + +/* From PBRT: core/montecarlo.h */ +inline float VanDerCorput(uint32_t n, uint32_t scramble) +{ + n = (n << 16) | (n >> 16); + n = ((n & 0x00ff00ff) << 8) | ((n & 0xff00ff00) >> 8); + n = ((n & 0x0f0f0f0f) << 4) | ((n & 0xf0f0f0f0) >> 4); + n = ((n & 0x33333333) << 2) | ((n & 0xcccccccc) >> 2); + n = ((n & 0x55555555) << 1) | ((n & 0xaaaaaaaa) >> 1); + n ^= scramble; + return ((n >> 8) & 0xffffff) / float(1 << 24); +} +inline float Sobol2(uint32_t n, uint32_t scramble) +{ + for (uint32_t v = 1 << 31; n != 0; n >>= 1, v ^= v >> 1) + if (n & 0x1) + scramble ^= v; + return ((scramble >> 8) & 0xffffff) / float(1 << 24); +} + +static float precompute_sheen_E(float rough, float mu, float u1, float u2) +{ + PrincipledSheenBsdf bsdf; + bsdf.weight = one_float3(); + bsdf.type = CLOSURE_BSDF_PRINCIPLED_SHEEN_V2_ID; + bsdf.sample_weight = 1.0f; + bsdf.N = make_float3(0.0f, 0.0f, 1.0f); + bsdf.roughness = sqr(rough); + + float3 eval, omega_in, domega_in_dx, domega_in_dy; + float pdf = 0.0f; + bsdf_principled_sheen_sample((ShaderClosure *)&bsdf, + make_float3(0.0f, 0.0f, 1.0f), + make_float3(sqrtf(1.0f - sqr(mu)), 0.0f, mu), + zero_float3(), + zero_float3(), + u1, + u2, + &eval, + &omega_in, + &domega_in_dx, + &domega_in_dy, + &pdf); + if (pdf != 0.0f) { + return clamp(average(eval) / pdf, 0.0f, 1e5f); + } + return 0.0f; +} + +static float precompute_clearcoat_E(float rough, float mu, float u1, float u2) +{ + MicrofacetBsdf bsdf; + bsdf.weight = one_float3(); + bsdf.type = CLOSURE_BSDF_MICROFACET_GGX_CLEARCOAT_V2_ID; + bsdf.sample_weight = 1.0f; + bsdf.N = make_float3(0.0f, 0.0f, 1.0f); + bsdf.alpha_x = bsdf.alpha_y = sqr(rough); + bsdf.ior = 1.5f; + bsdf.extra = nullptr; + bsdf.T = make_float3(1.0f, 0.0f, 0.0f); + + /* Account for the albedo scaling that the closure performs. + * Dependency warning - this relies on the ggx_E and ggx_E_avg lookup tables! */ + float E = microfacet_ggx_E(mu, rough), E_avg = microfacet_ggx_E_avg(rough); + float Fss = dielectric_fresnel_Fss(1.5f); + float Fms = Fss * E_avg / (1.0f - Fss * (1.0f - E_avg)); + float albedo_scale = 1.0f + Fms * ((1.0f - E) / E); + + float3 eval, omega_in, domega_in_dx, domega_in_dy; + float pdf = 0.0f; + bsdf_microfacet_ggx_sample((ShaderClosure *)&bsdf, + make_float3(0.0f, 0.0f, 1.0f), + make_float3(sqrtf(1.0f - sqr(mu)), 0.0f, mu), + zero_float3(), + zero_float3(), + u1, + u2, + &eval, + &omega_in, + &domega_in_dx, + &domega_in_dy, + &pdf); + if (pdf != 0.0f) { + /* Encode relative to macrosurface Fresnel, saves resolution. + * TODO: Worth the extra evaluation? */ + return albedo_scale * (average(eval) / pdf) / fresnel_dielectric_cos(mu, 1.5f); + } + return 0.0f; +} + +static float precompute_ggx_E(float rough, float mu, float u1, float u2) +{ + MicrofacetBsdf bsdf; + bsdf.weight = one_float3(); + bsdf.type = CLOSURE_BSDF_MICROFACET_GGX_ID; + bsdf.sample_weight = 1.0f; + bsdf.N = make_float3(0.0f, 0.0f, 1.0f); + bsdf.alpha_x = bsdf.alpha_y = sqr(rough); + bsdf.ior = 1.0f; + bsdf.extra = nullptr; + bsdf.T = make_float3(1.0f, 0.0f, 0.0f); + + float3 eval, omega_in, domega_in_dx, domega_in_dy; + float pdf = 0.0f; + bsdf_microfacet_ggx_sample((ShaderClosure *)&bsdf, + make_float3(0.0f, 0.0f, 1.0f), + make_float3(sqrtf(1.0f - sqr(mu)), 0.0f, mu), + zero_float3(), + zero_float3(), + u1, + u2, + &eval, + &omega_in, + &domega_in_dx, + &domega_in_dy, + &pdf); + if (pdf != 0.0f) { + return average(eval) / pdf; + } + return 0.0f; +} + +static float precompute_ggx_refract_E(float rough, float mu, float eta, float u1, float u2) +{ + MicrofacetBsdf bsdf; + bsdf.weight = one_float3(); + bsdf.type = CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID; + bsdf.sample_weight = 1.0f; + bsdf.N = make_float3(0.0f, 0.0f, 1.0f); + bsdf.alpha_x = bsdf.alpha_y = sqr(rough); + bsdf.ior = eta; + bsdf.extra = nullptr; + bsdf.T = make_float3(1.0f, 0.0f, 0.0f); + + float3 eval, omega_in, domega_in_dx, domega_in_dy; + float pdf = 0.0f; + bsdf_microfacet_ggx_sample((ShaderClosure *)&bsdf, + make_float3(0.0f, 0.0f, 1.0f), + make_float3(sqrtf(1.0f - sqr(mu)), 0.0f, mu), + zero_float3(), + zero_float3(), + u1, + u2, + &eval, + &omega_in, + &domega_in_dx, + &domega_in_dy, + &pdf); + if (pdf != 0.0f) { + return average(eval) / pdf; + } + return 0.0f; +} + +static float precompute_ggx_glass_E(float rough, float mu, float eta, float u1, float u2) +{ + MicrofacetBsdf bsdf; + bsdf.weight = one_float3(); + bsdf.type = CLOSURE_BSDF_MICROFACET_MULTI_GGX_GLASS_ID; + bsdf.sample_weight = 1.0f; + bsdf.N = make_float3(0.0f, 0.0f, 1.0f); + bsdf.alpha_x = bsdf.alpha_y = sqr(rough); + bsdf.ior = eta; + bsdf.extra = nullptr; + bsdf.T = make_float3(1.0f, 0.0f, 0.0f); + + float3 eval, omega_in, domega_in_dx, domega_in_dy; + float pdf = 0.0f; + bsdf_microfacet_ggx_glass_sample((ShaderClosure *)&bsdf, + make_float3(0.0f, 0.0f, 1.0f), + make_float3(sqrtf(1.0f - sqr(mu)), 0.0f, mu), + zero_float3(), + zero_float3(), + u1, + u2, + &eval, + &omega_in, + &domega_in_dx, + &domega_in_dy, + &pdf); + if (pdf != 0.0f) { + return average(eval) / pdf; + } + return 0.0f; +} + +static float precompute_ggx_dielectric_E(float rough, float mu, float eta, float u1, float u2) +{ + { + /* Reparametrize based on macrosurface fresnel to get more resolution into areas where + * the Fresnel curve is rapidly changing. Particularly important for eta<1 due to the TIR edge. + * However, in the eta<1 case, the entire TIR area would be compressed down to a point, which + * is an issue since there are changes in that range at higher roughnesses. + * Therefore, the remapping is blended with the identity function for a compromise. + */ + float F0 = fresnel_dielectric_cos(1.0f, eta); + auto get_remap = [eta, F0](float x) { + return mix(x, inverse_lerp(1.0f, F0, fresnel_dielectric_cos(x, eta)), 0.5f); + }; + + float remap_target = mu; + float start = 0.0f, end = 1.0f; + while (end - start > 1e-7f) { + mu = (end + start) * 0.5f; + if (get_remap(mu) > remap_target) { + end = mu; + } + else { + start = mu; + } + } + } + + MicrofacetExtrav2 extra; + MicrofacetBsdf bsdf; + bsdf.weight = one_float3(); + bsdf.type = CLOSURE_BSDF_MICROFACET_GGX_FRESNEL_V2_ID; + bsdf.sample_weight = 1.0f; + bsdf.N = make_float3(0.0f, 0.0f, 1.0f); + bsdf.alpha_x = bsdf.alpha_y = sqr(rough); + bsdf.ior = eta; + bsdf.extra = (MicrofacetExtra *)&extra; + bsdf.T = make_float3(1.0f, 0.0f, 0.0f); + extra.metallic = extra.metal_base = extra.metal_edge_factor = zero_float3(); + + /* Dependency warning - this relies on the ggx_E and ggx_E_avg lookup tables! */ + float E = microfacet_ggx_E(mu, rough), E_avg = microfacet_ggx_E_avg(rough); + float Fss = dielectric_fresnel_Fss(eta); + float Fms = Fss * E_avg / (1.0f - Fss * (1.0f - E_avg)); + extra.dielectric = 1.0f + Fms * ((1.0f - E) / E); + + float3 eval, omega_in, domega_in_dx, domega_in_dy; + float pdf = 0.0f; + bsdf_microfacet_ggx_sample((ShaderClosure *)&bsdf, + make_float3(0.0f, 0.0f, 1.0f), + make_float3(sqrtf(1.0f - sqr(mu)), 0.0f, mu), + zero_float3(), + zero_float3(), + u1, + u2, + &eval, + &omega_in, + &domega_in_dx, + &domega_in_dy, + &pdf); + if (pdf != 0.0f) { + return average(eval) / pdf; + } + return 0.0f; +} + +struct PrecomputeTerm { + int dim, samples, res; + std::function<float(float, float, float, float, float)> evaluation; +}; + +bool cycles_precompute(std::string name); +bool cycles_precompute(std::string name) +{ + std::map<string, PrecomputeTerm> precompute_terms; + precompute_terms["sheen_E"] = { + 2, 1 << 23, 32, [](float rough, float mu, float ior, float u1, float u2) { + return precompute_sheen_E(rough, mu, u1, u2); + }}; + precompute_terms["clearcoat_E"] = { + 2, 1 << 23, 16, [](float rough, float mu, float ior, float u1, float u2) { + return precompute_clearcoat_E(rough, mu, u1, u2); + }}; + precompute_terms["ggx_E"] = { + 2, 1 << 23, 32, [](float rough, float mu, float ior, float u1, float u2) { + return precompute_ggx_E(rough, mu, u1, u2); + }}; + precompute_terms["ggx_E_avg"] = { + 1, 1 << 23, 32, [](float rough, float mu, float ior, float u1, float u2) { + return 2.0f * mu * precompute_ggx_E(rough, mu, u1, u2); + }}; + precompute_terms["ggx_glass_E"] = { + 3, 1 << 20, 16, [](float rough, float mu, float ior, float u1, float u2) { + return precompute_ggx_glass_E(rough, mu, ior, u1, u2); + }}; + precompute_terms["ggx_glass_inv_E"] = { + 3, 1 << 20, 16, [](float rough, float mu, float ior, float u1, float u2) { + return precompute_ggx_glass_E(rough, mu, 1.0f / ior, u1, u2); + }}; + precompute_terms["ggx_refract_E"] = { + 3, 1 << 20, 16, [](float rough, float mu, float ior, float u1, float u2) { + return precompute_ggx_refract_E(rough, mu, ior, u1, u2); + }}; + precompute_terms["ggx_refract_inv_E"] = { + 3, 1 << 20, 16, [](float rough, float mu, float ior, float u1, float u2) { + return precompute_ggx_refract_E(rough, mu, 1.0f / ior, u1, u2); + }}; + precompute_terms["ggx_dielectric_E"] = { + 3, 1 << 20, 16, [](float rough, float mu, float ior, float u1, float u2) { + return precompute_ggx_dielectric_E(rough, mu, ior, u1, u2); + }}; + // TODO: Consider more X resolution for this table. + precompute_terms["ggx_dielectric_inv_E"] = { + 3, 1 << 20, 16, [](float rough, float mu, float ior, float u1, float u2) { + return precompute_ggx_dielectric_E(rough, mu, 1.0f / ior, u1, u2); + }}; + + if (precompute_terms.count(name) == 0) { + return false; + } + + const PrecomputeTerm &term = precompute_terms[name]; + + const int samples = term.samples; + const int res = term.res; + const int nz = (term.dim > 2) ? res : 1, ny = res, nx = (term.dim > 1) ? res : 1; + + if (nz > 1) { + std::cout << "static const float table_" << name << "[" << nz << "][" << ny << "][" << nx + << "] = {" << std::endl; + } + for (int z = 0; z < nz; z++) { + float *data = new float[nx * ny]; + parallel_for(0, ny, [&](int64_t y) { + for (int x = 0; x < nx; x++) { + uint rng = hash_uint2(x, y); + uint scramble1 = lcg_step_uint(&rng), scramble2 = lcg_step_uint(&rng); + double sum = 0.0; + for (int i = 0; i < samples; i++) { + float rough = 1.0f - (float(y) + lcg_step_float(&rng)) / float(ny); + float mu = (float(x) + lcg_step_float(&rng)) / float(nx); + float ior = (float(z) + lcg_step_float(&rng)) / float(nz); + /* Encode IOR remapped as sqrt(0.5*(IOR-1)) for more resolution at the start, where most + * of the changes happen (also places the most common range around 1.5 in the center) */ + ior = 1.0f + 2.0f * sqr(ior); + float u1 = VanDerCorput(i, scramble1); + float u2 = Sobol2(i, scramble2); + + float value = term.evaluation(rough, mu, ior, u1, u2); + if (isnan(value)) { + value = 0.0f; + } + sum += (double)value; + } + data[y * nx + x] = float(sum / double(samples)); + } + }); + + string filename = name; + if (nz > 1) { + filename += string_printf("_%02d", z); + std::cout << " {" << std::endl; + } + else { + std::cout << "static const float table_" << name << "[" << ny << "][" << nx << "] = {" + << std::endl; + } + + for (int y = 0; y < ny; y++) { + std::cout << " {"; + for (int x = 0; x < nx; x++) { + std::cout << data[y * nx + x] << ((x + 1 == nx) ? "f" : "f, "); + } + std::cout << ((y + 1 == ny) ? "}" : "},") << std::endl; + } + if (nz > 1) { + std::cout << ((z + 1 == nz) ? " }" : " },") << std::endl; + } + else { + std::cout << "};" << std::endl; + } + + FILE *f = fopen((filename + ".pfm").c_str(), "w"); + fprintf(f, "Pf\n%d %d\n-1.0\n", nx, ny); + fwrite(data, sizeof(float), nx * ny, f); + fclose(f); + } + if (nz > 1) { + std::cout << "};" << std::endl; + } + + return true; +} + +CCL_NAMESPACE_END |