yocto_extension.cpp

//
// Implementation for Yocto/Extension.
//

//
// LICENSE:
//
// Copyright (c) 2020 -- 2020 Fabio Pellacini
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
//

#include "yocto_extension.h"

#include <atomic>
#include <deque>
#include <future>
#include <memory>
#include <mutex>
#include <stdexcept>

using namespace std::string_literals;

// -----------------------------------------------------------------------------
// MATH FUNCTIONS
// -----------------------------------------------------------------------------
namespace yocto::extension {

// import math symbols for use
using math::abs;
using math::acos;
using math::atan2;
using math::clamp;
using math::cos;
using math::exp;
using math::flt_eps;
using math::flt_max;
using math::fmod;
using math::fresnel_conductor;
using math::fresnel_dielectric;
using math::identity3x3f;
using math::invalidb3f;
using math::log;
using math::luminance;
using math::make_rng;
using math::max;
using math::min;
using math::pif;
using math::pow;
using math::rand1f;
using math::rand2f;
using math::sample_discrete_cdf;
using math::sample_discrete_cdf_pdf;
using math::sample_sphere;
using math::sample_sphere_pdf;
using math::sample_uniform;
using math::sample_uniform_pdf;
using math::sin;
using math::sqrt;
using math::zero2f;
using math::zero2i;
using math::zero3f;
using math::zero3i;
using math::zero4f;
using math::zero4i;

}  // namespace yocto::extension

// -----------------------------------------------------------------------------
// IMPLEMENTATION FOR EXTENSION
// -----------------------------------------------------------------------------
namespace yocto::extension {

inline const float sqrt_pi_over_8f = 0.626657069f;

inline float sqr(float v) { return v * v; }
inline vec3f sqr(vec3f v) { return v * v; }

template <int n>
static float pow(float v) {
  auto n2 = pow<n / 2>(v);
  return n2 * n2 * pow<n & 1>(v);
}

template <>
inline float pow<1>(float v) {
  return v;
}

template <>
inline float pow<0>(float v) {
  return 1;
}

inline float safe_asin(float x) { return asin(clamp(x, -1.0f, 1.0f)); }

inline float safe_sqrt(float x) { return sqrt(max(0.0f, x)); }

static vec3f sigma_a_from_concentration(float ce, float cp) {
  auto eumelanin_sigma_a   = vec3f{0.419f, 0.697f, 1.37f};
  auto pheomelanin_sigma_a = vec3f{0.187f, 0.4f, 1.05f};
  return ce * eumelanin_sigma_a + cp * pheomelanin_sigma_a;
}

static vec3f sigma_a_from_reflectance(const vec3f& c, float beta_n) {
  return sqr(log(c) / (5.969f - 0.215f * beta_n + 2.532f * sqr(beta_n) -
                          10.73f * pow<3>(beta_n) + 5.574f * pow<4>(beta_n) +
                          0.245f * pow<5>(beta_n)));
}

hair_brdf eval_hair_brdf(const hair_material& material, float v,
    const vec3f& normal, const vec3f& tangent) {
  auto brdf = hair_brdf{};

  if (material.sigma_a) {
    brdf.sigma_a = material.sigma_a;
  } else if (material.color) {
    brdf.sigma_a = sigma_a_from_reflectance(material.color, material.beta_n);
  } else if (material.eumelanin || material.pheomelanin) {
    brdf.sigma_a = sigma_a_from_concentration(
        material.eumelanin, material.pheomelanin);
  }

  auto beta_m = material.beta_m;
  auto beta_n = material.beta_n;
  brdf.alpha  = material.alpha;
  brdf.eta    = material.eta;

#ifdef YOCTO_EMBREE
  brdf.h = v;
#else
  brdf.h = -1 + 2 * v;
#endif

  brdf.gamma_o = safe_asin(brdf.h);

  // Compute longitudinal variance from $\beta_m$
  brdf.v[0] = sqr(
      0.726f * beta_m + 0.812f * sqr(beta_m) + 3.7f * pow<20>(beta_m));
  brdf.v[1] = 0.25f * brdf.v[0];
  brdf.v[2] = 4 * brdf.v[0];
  for (auto p = 3; p <= p_max; p++) brdf.v[p] = brdf.v[2];

  // Compute azimuthal logistic scale factor from $\beta_n$
  brdf.s = sqrt_pi_over_8f *
           (0.265f * beta_n + 1.194f * sqr(beta_n) + 5.372f * pow<22>(beta_n));

  // Compute $\alpha$ terms for hair scales
  brdf.sin_2k_alpha[0] = sin(pif / 180 * brdf.alpha);
  brdf.cos_2k_alpha[0] = safe_sqrt(1 - sqr(brdf.sin_2k_alpha[0]));
  for (auto i = 1; i < 3; i++) {
    brdf.sin_2k_alpha[i] = 2 * brdf.cos_2k_alpha[i - 1] *
                           brdf.sin_2k_alpha[i - 1];
    brdf.cos_2k_alpha[i] = sqr(brdf.cos_2k_alpha[i - 1]) -
                           sqr(brdf.sin_2k_alpha[i - 1]);
  }

  brdf.world_to_brdf = inverse(frame_fromzx(zero3f, normal, tangent));

  return brdf;
}

inline float i0(float x) {
  float   val   = 0;
  float   x2i   = 1;
  int64_t ifact = 1;
  int     i4    = 1;
  // I0(x) \approx Sum_i x^(2i) / (4^i (i!)^2)
  for (int i = 0; i < 10; i++) {
    if (i > 1) ifact *= i;
    val += x2i / (i4 * ifact * ifact);
    x2i *= x * x;
    i4 *= 4;
  }
  return val;
}

inline float log_i0(float x) {
  if (x > 12)
    return x + 0.5f * (-log(2 * pif) + log(1 / x) + 1 / (8 * x));
  else
    return log(i0(x));
}

static float mp(float cos_theta_i, float cos_theta_o, float sin_theta_i,
    float sin_theta_o, float v) {
  auto a = cos_theta_i * cos_theta_o / v;
  auto b = sin_theta_i * sin_theta_o / v;
  return (v <= 0.1f) ? (exp(log_i0(a) - b - 1 / v + 0.6931f + log(1 / (2 * v))))
                     : (exp(-b) * i0(a)) / (sinh(1 / v) * 2 * v);
}

static std::array<vec3f, p_max + 1> ap(
    float cos_theta_o, float eta, float h, const vec3f& T) {
  auto ap = std::array<vec3f, p_max + 1>{};
  // Compute $p=0$ attenuation at initial cylinder intersection
  auto cos_gamma_o = safe_sqrt(1 - h * h);
  auto cos_theta   = cos_theta_o * cos_gamma_o;

  // fresnel_dielectric(eta, cos_theta): we hack function's parameters bulding
  // two vectors s.t. their dot product gives exactly cos_theta
  auto f = fresnel_dielectric(eta, {0, 0, 1}, {0, 0, cos_theta});
  ap[0]  = vec3f(f);

  // Compute $p=1$ attenuation term
  ap[1] = sqr(1 - f) * T;

  // Compute attenuation terms up to $p=_pMax_$
  for (auto p = 2; p < p_max; p++) ap[p] = ap[p - 1] * T * f;

  // Compute attenuation term accounting for remaining orders of scattering
  ap[p_max] = ap[p_max - 1] * f * T / (vec3f(1.f) - T * f);
  return ap;
}

inline float phi(int p, float gamma_o, float gamma_t) {
  return 2 * p * gamma_t - 2 * gamma_o + p * pif;
}

inline float logistic(float x, float s) {
  x = abs(x);
  return exp(-x / s) / (s * sqr(1 + exp(-x / s)));
}

inline float logistic_cdf(float x, float s) { return 1 / (1 + exp(-x / s)); }

inline float trimmed_logistic(float x, float s, float a, float b) {
  return logistic(x, s) / (logistic_cdf(b, s) - logistic_cdf(a, s));
}

inline float np(float phi, int p, float s, float gamma_o, float gamma_t) {
  auto dphi = phi - ext::phi(p, gamma_o, gamma_t);
  // Remap _dphi_ to $[-\pi,\pi]$
  while (dphi > pif) dphi -= 2 * pif;
  while (dphi < -pif) dphi += 2 * pif;
  return trimmed_logistic(dphi, s, -pif, pif);
}

vec3f eval_hair_scattering(
    const hair_brdf& brdf, const vec3f& outgoing_, const vec3f& incoming_) {
  auto sigma_a       = brdf.sigma_a;
  auto eta           = brdf.eta;
  auto h             = brdf.h;
  auto gamma_o       = brdf.gamma_o;
  auto v             = brdf.v;
  auto s             = brdf.s;
  auto sin_2k_alpha  = brdf.sin_2k_alpha;
  auto cos_2k_alpha  = brdf.cos_2k_alpha;
  auto world_to_brdf = brdf.world_to_brdf;

  auto outgoing = transform_direction(world_to_brdf, outgoing_);
  auto incoming = transform_direction(world_to_brdf, incoming_);

  // Compute hair coordinate system terms related to _wo_
  auto sin_theta_o = outgoing.x;
  auto cos_theta_o = safe_sqrt(1 - sqr(sin_theta_o));
  auto phi_o       = atan2(outgoing.z, outgoing.y);

  // Compute hair coordinate system terms related to _wi_
  auto sin_theta_i = incoming.x;
  auto cos_theta_i = safe_sqrt(1 - sqr(sin_theta_i));
  auto phi_i       = atan2(incoming.z, incoming.y);

  // Compute $\cos \thetat$ for refracted ray
  auto sin_theta_t = sin_theta_o / eta;
  auto cos_theta_t = safe_sqrt(1 - sqr(sin_theta_t));

  // Compute $\gammat$ for refracted ray
  auto etap        = sqrt(eta * eta - sqr(sin_theta_o)) / cos_theta_o;
  auto sin_gamma_t = h / etap;
  auto cos_gamma_t = safe_sqrt(1 - sqr(sin_gamma_t));
  auto gamma_t     = safe_asin(sin_gamma_t);

  // Compute the transmittance _T_ of a single path through the cylinder
  auto T = exp(-sigma_a * (2 * cos_gamma_t / cos_theta_t));

  // Evaluate hair BSDF
  auto phi  = phi_i - phi_o;
  auto ap   = ext::ap(cos_theta_o, eta, h, T);
  auto fsum = zero3f;
  for (auto p = 0; p < p_max; p++) {
    // Compute $\sin \thetao$ and $\cos \thetao$ terms accounting for scales
    auto sin_theta_op = 0.0f;
    auto cos_theta_op = 0.0f;
    if (p == 0) {
      sin_theta_op = sin_theta_o * cos_2k_alpha[1] -
                     cos_theta_o * sin_2k_alpha[1];
      cos_theta_op = cos_theta_o * cos_2k_alpha[1] +
                     sin_theta_o * sin_2k_alpha[1];
    }

    // Handle remainder of $p$ values for hair scale tilt
    else if (p == 1) {
      sin_theta_op = sin_theta_o * cos_2k_alpha[0] +
                     cos_theta_o * sin_2k_alpha[0];
      cos_theta_op = cos_theta_o * cos_2k_alpha[0] -
                     sin_theta_o * sin_2k_alpha[0];
    } else if (p == 2) {
      sin_theta_op = sin_theta_o * cos_2k_alpha[2] +
                     cos_theta_o * sin_2k_alpha[2];
      cos_theta_op = cos_theta_o * cos_2k_alpha[2] -
                     sin_theta_o * sin_2k_alpha[2];
    } else {
      sin_theta_op = sin_theta_o;
      cos_theta_op = cos_theta_o;
    }

    // Handle out-of-range $\cos \thetao$ from scale adjustment
    cos_theta_op = abs(cos_theta_op);
    fsum += mp(cos_theta_i, cos_theta_op, sin_theta_i, sin_theta_op, v[p]) *
            ap[p] * np(phi, p, s, gamma_o, gamma_t);
  }

  // Compute contribution of remaining terms after _pMax_
  fsum += mp(cos_theta_i, cos_theta_o, sin_theta_i, sin_theta_o, v[p_max]) *
          ap[p_max] / (2 * pif);
  // if (abs(incoming.z) > 0)
  //   fsum /= abs(incoming.z);
  return fsum;
}

// https://fgiesen.wordpress.com/2009/12/13/decoding-morton-codes/
static uint32_t compact1by1(uint32_t x) {
  // x = -f-e -d-c -b-a -9-8 -7-6 -5-4 -3-2 -1-0
  x &= 0x55555555;
  // x = --fe --dc --ba --98 --76 --54 --32 --10
  x = (x ^ (x >> 1)) & 0x33333333;
  // x = ---- fedc ---- ba98 ---- 7654 ---- 3210
  x = (x ^ (x >> 2)) & 0x0f0f0f0f;
  // x = ---- ---- fedc ba98 ---- ---- 7654 3210
  x = (x ^ (x >> 4)) & 0x00ff00ff;
  // x = ---- ---- ---- ---- fedc ba98 7654 3210
  x = (x ^ (x >> 8)) & 0x0000ffff;
  return x;
}

static vec2f demux_float(float f) {
  uint64_t v       = f * (1ull << 32);
  uint32_t bits[2] = {compact1by1(v), compact1by1(v >> 1)};
  return {bits[0] / float(1 << 16), bits[1] / float(1 << 16)};
}

static float sample_trimmed_logistic(float u, float s, float a, float b) {
  auto k = logistic_cdf(b, s) - logistic_cdf(a, s);
  auto x = -s * log(1 / (u * k + logistic_cdf(a, s)) - 1);
  return clamp(x, a, b);
}

static std::array<float, p_max + 1> compute_ap_pdf(
    const hair_brdf& brdf, float cos_theta_o) {
  auto sigma_a = brdf.sigma_a;
  auto eta     = brdf.eta;
  auto h       = brdf.h;

  // Compute array of $A_p$ values for _cosThetaO_
  auto sin_theta_o = safe_sqrt(1 - cos_theta_o * cos_theta_o);

  // Compute $\cos \thetat$ for refracted ray
  auto sin_theta_t = sin_theta_o / eta;
  auto cos_theta_t = safe_sqrt(1 - sqr(sin_theta_t));

  // Compute $\gammat$ for refracted ray
  auto etap        = sqrt(eta * eta - sqr(sin_theta_o)) / cos_theta_o;
  auto sin_gamma_t = h / etap;
  auto cos_gamma_t = safe_sqrt(1 - sqr(sin_gamma_t));

  // Compute the transmittance _T_ of a single path through the cylinder
  auto T  = exp(-sigma_a * (2 * cos_gamma_t / cos_theta_t));
  auto ap = ext::ap(cos_theta_o, eta, h, T);

  // Compute $A_p$ PDF from individual $A_p$ terms
  auto ap_pdf = std::array<float, p_max + 1>{};
  auto sum_y  = 0.0f;
  for (auto i = 0; i <= p_max; i++) {
    sum_y += luminance(ap[i]);
  }
  for (auto i = 0; i <= p_max; i++) {
    ap_pdf[i] = luminance(ap[i]) / sum_y;
  }
  return ap_pdf;
}

vec3f sample_hair_scattering(
    const hair_brdf& brdf, const vec3f& outgoing_, const vec2f& rn) {
  auto eta           = brdf.eta;
  auto h             = brdf.h;
  auto gamma_o       = brdf.gamma_o;
  auto v             = brdf.v;
  auto s             = brdf.s;
  auto sin_2k_alpha  = brdf.sin_2k_alpha;
  auto cos_2k_alpha  = brdf.cos_2k_alpha;
  auto world_to_brdf = brdf.world_to_brdf;

  auto outgoing = transform_direction(world_to_brdf, outgoing_);

  // Compute hair coordinate system terms related to _wo_
  auto sin_theta_o = outgoing.x;
  auto cos_theta_o = safe_sqrt(1 - sqr(sin_theta_o));
  auto phi_o       = atan2(outgoing.z, outgoing.y);

  // Derive four random samples from _u2_
  auto u = std::array<vec2f, 2>{demux_float(rn.x), demux_float(rn.y)};

  // Determine which term $p$ to sample for hair scattering
  auto ap_pdf = compute_ap_pdf(brdf, cos_theta_o);
  auto p      = 0;
  for (p = 0; p < p_max; p++) {
    if (u[0][0] < ap_pdf[p]) break;
    u[0][0] -= ap_pdf[p];
  }

  // Rotate $\sin \thetao$ and $\cos \thetao$ to account for hair scale tilt
  auto sin_theta_op = 0.0f;
  auto cos_theta_op = 0.0f;
  if (p == 0) {
    sin_theta_op = sin_theta_o * cos_2k_alpha[1] -
                   cos_theta_o * sin_2k_alpha[1];
    cos_theta_op = cos_theta_o * cos_2k_alpha[1] +
                   sin_theta_o * sin_2k_alpha[1];
  } else if (p == 1) {
    sin_theta_op = sin_theta_o * cos_2k_alpha[0] +
                   cos_theta_o * sin_2k_alpha[0];
    cos_theta_op = cos_theta_o * cos_2k_alpha[0] -
                   sin_theta_o * sin_2k_alpha[0];
  } else if (p == 2) {
    sin_theta_op = sin_theta_o * cos_2k_alpha[2] +
                   cos_theta_o * sin_2k_alpha[2];
    cos_theta_op = cos_theta_o * cos_2k_alpha[2] -
                   sin_theta_o * sin_2k_alpha[2];
  } else {
    sin_theta_op = sin_theta_o;
    cos_theta_op = cos_theta_o;
  }

  // Sample $M_p$ to compute $\thetai$
  u[1][0]          = max(u[1][0], 1e-5f);
  auto cos_theta   = 1 + v[p] * log(u[1][0] + (1 - u[1][0]) * exp(-2 / v[p]));
  auto sin_theta   = safe_sqrt(1 - sqr(cos_theta));
  auto cos_phi     = cos(2 * pif * u[1][1]);
  auto sin_theta_i = -cos_theta * sin_theta_op +
                     sin_theta * cos_phi * cos_theta_op;
  auto cos_theta_i = safe_sqrt(1 - sqr(sin_theta_i));

  // Sample $N_p$ to compute $\Delta\phi$

  // Compute $\gammat$ for refracted ray
  auto etap        = sqrt(eta * eta - sqr(sin_theta_o)) / cos_theta_o;
  auto sin_gamma_t = h / etap;
  auto gamma_t     = safe_asin(sin_gamma_t);
  auto dphi        = 0.0f;
  if (p < p_max)
    dphi = phi(p, gamma_o, gamma_t) +
           sample_trimmed_logistic(u[0][1], s, -pif, pif);
  else
    dphi = 2 * pif * u[0][1];

  // Compute _wi_ from sampled hair scattering angles
  auto phi_i = phi_o + dphi;

  auto incoming = vec3f{
      sin_theta_i, cos_theta_i * cos(phi_i), cos_theta_i * sin(phi_i)};
  return transform_direction(inverse(world_to_brdf), incoming);
}

float sample_hair_scattering_pdf(
    const hair_brdf& brdf, const vec3f& outgoing_, const vec3f& incoming_) {
  auto eta           = brdf.eta;
  auto h             = brdf.h;
  auto gamma_o       = brdf.gamma_o;
  auto v             = brdf.v;
  auto s             = brdf.s;
  auto sin_2k_alpha  = brdf.sin_2k_alpha;
  auto cos_2k_alpha  = brdf.cos_2k_alpha;
  auto world_to_brdf = brdf.world_to_brdf;

  auto outgoing = transform_direction(world_to_brdf, outgoing_);
  auto incoming = transform_direction(world_to_brdf, incoming_);

  // Compute hair coordinate system terms related to _wo_
  auto sin_theta_o = outgoing.x;
  auto cos_theta_o = safe_sqrt(1 - sqr(sin_theta_o));
  auto phi_o       = atan2(outgoing.z, outgoing.y);

  // Compute hair coordinate system terms related to _wi_
  auto sin_theta_i = incoming.x;
  auto cos_theta_i = safe_sqrt(1 - sqr(sin_theta_i));
  auto phi_i       = atan2(incoming.z, incoming.y);

  // Compute $\gammat$ for refracted ray
  auto etap        = sqrt(eta * eta - sqr(sin_theta_o)) / cos_theta_o;
  auto sin_gamma_t = h / etap;
  auto gamma_t     = safe_asin(sin_gamma_t);

  // Compute PDF for $A_p$ terms
  auto ap_pdf = compute_ap_pdf(brdf, cos_theta_o);

  // Compute PDF sum for hair scattering events
  auto phi = phi_i - phi_o;
  auto pdf = 0.0f;
  for (auto p = 0; p < p_max; p++) {
    // Compute $\sin \thetao$ and $\cos \thetao$ terms accounting for scales
    auto sin_theta_op = 0.0f;
    auto cos_theta_op = 0.0f;
    if (p == 0) {
      sin_theta_op = sin_theta_o * cos_2k_alpha[1] -
                     cos_theta_o * sin_2k_alpha[1];
      cos_theta_op = cos_theta_o * cos_2k_alpha[1] +
                     sin_theta_o * sin_2k_alpha[1];
    }

    // Handle remainder of $p$ values for hair scale tilt
    else if (p == 1) {
      sin_theta_op = sin_theta_o * cos_2k_alpha[0] +
                     cos_theta_o * sin_2k_alpha[0];
      cos_theta_op = cos_theta_o * cos_2k_alpha[0] -
                     sin_theta_o * sin_2k_alpha[0];
    } else if (p == 2) {
      sin_theta_op = sin_theta_o * cos_2k_alpha[2] +
                     cos_theta_o * sin_2k_alpha[2];
      cos_theta_op = cos_theta_o * cos_2k_alpha[2] -
                     sin_theta_o * sin_2k_alpha[2];
    } else {
      sin_theta_op = sin_theta_o;
      cos_theta_op = cos_theta_o;
    }

    // Handle out-of-range $\cos \thetao$ from scale adjustment
    cos_theta_op = abs(cos_theta_op);
    pdf += mp(cos_theta_i, cos_theta_op, sin_theta_i, sin_theta_op, v[p]) *
           ap_pdf[p] * np(phi, p, s, gamma_o, gamma_t);
  }
  pdf += mp(cos_theta_i, cos_theta_o, sin_theta_i, sin_theta_o, v[p_max]) *
         ap_pdf[p_max] * (1 / (2 * pif));
  return pdf;
}

// HAIR SCATTERING TESTS

void white_furnace_test() {
  auto rng = make_rng(199382389514);
  auto wo  = sample_sphere(rand2f(rng));
  for (auto beta_m = 0.1f; beta_m < 1.0f; beta_m += 0.2f) {
    for (auto beta_n = 0.1f; beta_n < 1.0f; beta_n += 0.2f) {
      // Estimate reflected uniform incident radiance from hair
      auto sum   = zero3f;
      auto count = 300000;
      for (auto i = 0; i < count; i++) {
#ifdef YOCTO_EMBREE
        auto h = -1 + 2 * rand1f(rng);
#else
        auto h = rand1f(rng);
#endif
        // the original pbrt test fails when h == 0
        if (h == 0) h += flt_eps;

        auto mat   = hair_material{};
        mat.beta_m = beta_m;
        mat.beta_n = beta_n;
        mat.alpha  = 0;
        auto brdf  = eval_hair_brdf(mat, h, {0, 0, 1}, {1, 0, 0});
        auto wi    = sample_sphere(rand2f(rng));
        sum += eval_hair_scattering(brdf, wo, wi);
      }
      auto avg = luminance(sum) / (count * sample_sphere_pdf(wo));
      if (!(avg >= 0.95f && avg <= 1.05f))
        throw std::runtime_error("TEST FAILED!");
    }
  }
  printf("OK!\n");
  fflush(stdout);
}

void white_furnace_sampled_test() {
  auto rng = make_rng(199382389514);
  auto wo  = sample_sphere(rand2f(rng));
  for (auto beta_m = 0.1f; beta_m < 1.0f; beta_m += 0.2f) {
    for (auto beta_n = 0.1f; beta_n < 1.0f; beta_n += 0.2f) {
      auto sum   = zero3f;
      auto count = 300000;
      for (auto i = 0; i < count; i++) {
#ifdef YOCTO_EMBREE
        auto h = -1 + 2 * rand1f(rng);
#else
        auto h = rand1f(rng);
#endif
        auto mat   = hair_material{};
        mat.beta_m = beta_m;
        mat.beta_n = beta_n;
        mat.alpha  = 0;
        auto brdf  = eval_hair_brdf(mat, h, {0, 0, 1}, {1, 0, 0});
        auto wi    = sample_hair_scattering(brdf, wo, rand2f(rng));
        auto f     = eval_hair_scattering(brdf, wo, wi);
        auto pdf   = sample_hair_scattering_pdf(brdf, wo, wi);
        if (pdf > 0) sum += f / pdf;
      }
      auto avg = luminance(sum) / (count);
      if (!(avg >= 0.99f && avg <= 1.01f))
        throw std::runtime_error("TEST FAILED!");
    }
  }
  printf("OK!\n");
  fflush(stdout);
}

void sampling_weights_test() {
  auto rng = make_rng(199382389514);
  for (auto beta_m = 0.1f; beta_m < 1.0f; beta_m += 0.2f) {
    for (auto beta_n = 0.4f; beta_n < 1.0f; beta_n += 0.2f) {
      auto count = 10000;
      for (auto i = 0; i < count; i++) {
        // Check _HairBSDF::Sample\_f()_ sample weight
#ifdef YOCTO_EMBREE
        auto h = -1 + 2 * rand1f(rng);
#else
        auto h = rand1f(rng);
#endif
        auto mat   = hair_material{};
        mat.beta_m = beta_m;
        mat.beta_n = beta_n;
        mat.alpha  = 0;
        auto brdf  = eval_hair_brdf(mat, h, {0, 0, 1}, {1, 0, 0});
        auto wo    = sample_sphere(rand2f(rng));
        auto wi    = sample_hair_scattering(brdf, wo, rand2f(rng));
        auto f     = eval_hair_scattering(brdf, wo, wi);
        auto pdf   = sample_hair_scattering_pdf(brdf, wo, wi);
        if (pdf > 0) {
          // Verify that hair BSDF sample weight is close to 1 for _wi_
          if (!(luminance(f) / pdf >= 0.999f && luminance(f) / pdf <= 1.001f))
            throw std::runtime_error("TEST FAILED!");
        }
      }
    }
  }
  printf("OK!\n");
  fflush(stdout);
}

void sampling_consistency_test() {
  auto rng = make_rng(199382389514);
  for (auto beta_m = 0.2f; beta_m < 1.0f; beta_m += 0.2f)
    for (auto beta_n = 0.4f; beta_n < 1.0f; beta_n += 0.2f) {
      // Declare variables for hair sampling test
      const auto count   = 64 * 1024;
      auto       sigma_a = vec3f(0.25f);
      auto       wo      = sample_sphere(rand2f(rng));
      auto       li = [](const vec3f& w) -> vec3f { return vec3f(w.z * w.z); };
      auto       f_importance = zero3f;
      auto       f_uniform    = zero3f;
      for (auto i = 0; i < count; i++) {
        // Compute estimates of scattered radiance for hair sampling test
#ifdef YOCTO_EMBREE
        auto h = -1 + 2 * rand1f(rng);
#else
        auto h = rand1f(rng);
#endif
        auto mat   = hair_material{};
        mat.beta_m = beta_m;
        mat.beta_n = beta_n;
        mat.alpha  = 0;
        auto brdf  = eval_hair_brdf(mat, h, {0, 0, 1}, {1, 0, 0});
        auto u     = rand2f(rng);
        auto wi    = sample_hair_scattering(brdf, wo, u);
        auto f     = eval_hair_scattering(brdf, wo, wi);
        auto pdf   = sample_hair_scattering_pdf(brdf, wo, wi);
        if (pdf > 0) f_importance += f * li(wi) / (count * pdf);
        wi = sample_sphere(u);
        f_uniform += eval_hair_scattering(brdf, wo, wi) * li(wi) /
                     (count * sample_sphere_pdf(wo));
      }
      // Verify consistency of estimated hair reflected radiance values
      auto err = abs(luminance(f_importance) - luminance(f_uniform)) /
                 luminance(f_uniform);
      if (err >= 0.05f) throw std::runtime_error("TEST FAILED!");
    }
  printf("OK!\n");
  fflush(stdout);
}

}  // namespace yocto::extension