#include "micrograin.h"
#include "sigma_tab.h"
//#pragma once

namespace LT_NAMESPACE {

    Float MicrograinMicrosurface::one_to_many(const Float& tau_0, const Float& sigma_) {
        Float rho = -std::log(1 - tau_0) / pi;
        return  std::exp(-rho * sigma_);
    }

    Float MicrograinMicrosurface::tau_v(const Float& tau_0, const vec3& wi_u) {
        return 1 - one_to_many(tau_0, sigma(wi_u));
    }

    Float MicrograinMicrosurface::D(const vec3& wh_u, const SurfaceInteraction& si) {
        Float tau = tau_0->eval(si);
        Float rho = -std::log(1 - tau) / pi;
        float D_ = rho * one_to_many(tau, sigma_base(wh_u)) / (tau);
        return D_;
    }

    Float MicrograinMicrosurface::D(const vec3& wh_u, const vec3& wi_u, const SurfaceInteraction& si) {
        Float tau = tau_0->eval(si);
        Float rho = -std::log(1 - tau) / pi;
        Float D_ = rho * one_to_many(tau, sigma_shadow(wh_u, wi_u) + sigma_base(wh_u)) * glm::max(glm::dot(wh_u, wi_u), 0.f) / (wi_u[2] * tau_v(tau,wi_u));;
        // Same as :
        // Float D_ = tau_0 * D(wh_u) * G1(wh_u,wi_u) * glm::max(glm::dot(wh_u,wi_u),0.f) / (wi_u[2] * tau_v(wi_u));
        assert(D_ >= 0);
        assert(D_ == D_);
        return D_;
    }

    Float MicrograinMicrosurface::pdf(const vec3& wh_u, const SurfaceInteraction& si)
    {
        return D(wh_u, si) * wh_u.z;
    }
    Float MicrograinMicrosurface::pdf(const vec3& wh_u, const vec3& wi_u, const SurfaceInteraction& si) {
        return D(wh_u, wi_u, si);
    }

    vec3 MicrograinMicrosurface::sample_D(Sampler& sampler, const SurfaceInteraction& si)
    {
        Float tau = tau_0->eval(si);
        float xi = sampler.next_float();
        // warp the uniform distribution based on the filling factor
        float Xi = std::log(1. - xi * tau) / std::log(1. - tau);

        return square_to_cosine_hemisphere(1-Xi, sampler.next_float());
    }


    
    inline Float f(const Float& phi, const Float& theta_i, const Float& theta_m, const Float phi_q) {
        Float num = std::max(std::sin(theta_i) * std::sin(theta_m) * std::cos(phi) + std::cos(theta_i) * std::cos(theta_m), 0.f);
        Float f_ = 0.5 * num / (std::sin(phi_q) * std::sin(theta_i) * std::sin(theta_m) + phi_q * std::cos(theta_i) * std::cos(theta_m));
        assert(f_ == f_);
        return f_;
    }


    inline Float f_prime(const Float& phi, const Float& theta_i, const Float& theta_m, const Float phi_q) {
        Float num = -std::sin(theta_i) * std::sin(theta_m) * std::sin(phi);
        return 0.5 * num / (std::sin(phi_q) * std::sin(theta_i) * std::sin(theta_m) + phi_q * std::cos(theta_i) * std::cos(theta_m));
    }

    inline Float F(const Float& phi, const Float& theta_i, const Float& theta_m, const Float phi_q) {
        Float num = std::sin(phi) * std::sin(theta_i) * std::sin(theta_m) + phi * std::cos(theta_i) * std::cos(theta_m);
        return 0.5 * (num / (std::sin(phi_q) * std::sin(theta_i) * std::sin(theta_m) + phi_q * std::cos(theta_i) * std::cos(theta_m)) + 1);
    }

    inline Float g(const Float& phi, const Float& theta_i, const Float& theta_m, const Float phi_q) {
        Float a = f(0, theta_i, theta_m, phi_q);
        Float b = f(phi_q, theta_i, theta_m, phi_q);
        Float res = ((b - a) / phi_q) * std::abs(phi) + a;
        return res / (phi_q * (a + b));
    }

    inline Float inv_G(const Float& u, const Float& theta_i, const Float& theta_m, const Float phi_q) {
        Float a = f(0, theta_i, theta_m, phi_q);
        Float b = f(phi_q, theta_i, theta_m, phi_q);
        if (std::abs(a - b) < 0.00001)
            return (2 * u - 1) * phi_q;

        Float Y = (2 * u - 1) * phi_q * (a + b);
        Float s = glm::sign(u - 0.5);
        Float x = -s * std::sqrt(std::max(a * a * phi_q * phi_q - s * a * phi_q * Y + s * b * phi_q * Y,0.f)) + s * a * phi_q;
        Float inv_G_ = x / (a - b);
        assert(inv_G_ == inv_G_);
        return inv_G_;
    }

    /*inline void sample_tent(const Float& u, const Float& a, const Float& b, const Float& phq, Float& x_a, Float& g_x_a) {
        if (std::abs(a - b) < 0.00001){
            x_a = (2 * u - 1) * phq;
            g_x_a = 1 / (2. * phq);
            return;
        }

        Float Y = (2 * u - 1) * phq * (a + b);
        Float s = glm::sign(u - 0.5);
        Float x = -s * std::sqrt(std::max(a * a * phq * phq - s * a * phq * Y + s * b * phq * Y, 0.f)) + s * a * phq;
        x_a = x / (a - b);
        g_x_a = (((b - a) / phq) * std::abs(x_a) + a) / (phq * (a + b));
        assert(x_a == x_a);
        assert(g_x_a == g_x_a);
    }*/

    void sample_tent(float xi, float phi_q, float h0, float& phi, float& pdf)
    {
        float h1 = 1. / phi_q - h0;
     
        if (std::abs(h0 - h1) < 0.00001) {
            phi = (2 * xi - 1) * phi_q;
            pdf = 1 / (2. * phi_q);
            return;
        }

        float cst = phi_q / (h1 - h0);
        float h1_sqr = h1 * h1;
        float h0_sqr = h0 * h0;
        float tmp = 0;
        if (xi < 0.5)
            tmp = std::sqrt(std::max(2. * h1_sqr * xi + h0_sqr * (1. - 2. * xi),0.)) - h1;
        else
            tmp = h1 - std::sqrt(std::max((h1 - h0) * (1. - 2. * xi) / phi_q + h1_sqr,0.));
        phi = tmp * cst;
        pdf = (h0 - h1) / phi_q * std::abs(phi) + h1;
        assert(phi == phi);
        assert(pdf == pdf);
    }


    inline Float tricut_remove_bias(const Float& u, const Float& v, const Float& x_a, const Float& f_x_a, const Float& F_x_a, const Float& fprime_x_a, const Float& g_x_a) {
        Float x_b = x_a + 2.0 * (u - F_x_a) / std::max(f_x_a, 0.000001f);

        Float w_a = f_x_a * f_x_a + 2.0 * (u - F_x_a) * fprime_x_a;
        Float w_b = 2.0 * f_x_a * g_x_a - w_a;
        if ((w_a < 0.000001f) && (w_b < 0.000001f))
            return x_a;

        Float t = v * (w_a + w_b) / (w_b + std::sqrt(std::max((1 - v) * w_b * w_b + v * w_a * w_a, 0.f)));

        Float x = x_a * t + x_b * (1 - t);
        assert(x == x);
        return x;
    }

    inline void sample_bias_phi(const Float& u, const Float& cos_thi, const Float& sin_thi, const Float& cos_th, const Float& sin_th, const Float& phq, Float& x_a, Float& f_x_a, Float& F_x_a, Float& fprime_x_a, Float& g_x_a) {
        #if 0
        x_a = inv_G(u, theta_i, theta_m, phq);
        f_x_a = f(x_a, theta_i, theta_m, phq);
        F_x_a = F(x_a, theta_i, theta_m, phq);
        fprime_x_a = f_prime(x_a, theta_i, theta_m, phq);
        g_x_a = g(x_a, theta_i, theta_m, phq);
        #else
        Float cos_phq = std::cos(phq);
        Float sin_phq = std::sin(phq);

        Float cos_thi_th = cos_thi * cos_th;
        Float sin_thi_th = sin_thi * sin_th;
        Float cst = 1./ (2. * (sin_phq * sin_thi_th + phq * cos_thi_th));
        Float h0 = cst * (cos_phq * sin_thi_th + cos_thi_th);
        sample_tent(u, phq, h0, x_a, g_x_a);
        
        Float sin_ph = sin(x_a), cos_ph = cos(x_a);
        F_x_a = 0.5 + cst * (sin_ph * sin_thi_th + x_a * cos_thi_th);
        f_x_a = cst * (cos_ph * sin_thi_th + cos_thi_th);
        fprime_x_a = -cst * sin_ph * sin_thi_th;

        #endif  
    }

    inline Float sample_phi_tricut(const Float& u, const Float& v, const Float& cos_thi, const Float& sin_thi, const Float& cos_th, const Float& sin_th, const Float phi_q) {    
        Float x_a, f_x_a, F_x_a, fprime_x_a, g_x_a;
        sample_bias_phi(u, cos_thi, sin_thi, cos_th, sin_th, phi_q, x_a, f_x_a, F_x_a, fprime_x_a, g_x_a);
        return tricut_remove_bias(u, v, x_a, f_x_a, F_x_a, fprime_x_a, g_x_a);
    }

    Float MicrograinMicrosurface::sample_phi(const Float& cos_th, const vec3& wi_u, Sampler& sampler, const SurfaceInteraction& si) {
        Float sin_thi = std::sqrt(1 - wi_u.z * wi_u.z);
        Float tan_thi = sin_thi / wi_u.z;
        Float sin_th  = std::sqrt(1 - cos_th * cos_th);
        Float tan_th  = sin_th  / cos_th;
        
        Float phq = pi;
        if(sin_thi > cos_th)
            phq = std::acos(-1. / (tan_thi * tan_th) );
        
        Float phi = sample_phi_tricut(sampler.next_float() , sampler.next_float(), wi_u.z, sin_thi, cos_th, sin_th, phq );
        
        assert(phi == phi);
        return phi;
    }

    Float MicrograinMicrosurface::sample_cos_theta(const vec3& wi_u, Sampler& sampler, const SurfaceInteraction& si) {
        Float tau = tau_0->eval(si);
        Float Xi = pi * std::log(1 - sampler.next_float() * tau_v(tau,wi_u)) / std::log(1 - tau);

        Float cos_theta = sigma_inv(Xi, wi_u);
    
        assert(cos_theta == cos_theta);
        return cos_theta;
    }

    vec3 MicrograinMicrosurface::sample_D(const vec3& wi_u, Sampler& sampler, const SurfaceInteraction& si)
    {
        Float cos_theta = sample_cos_theta(wi_u, sampler, si);
        Float sin_theta = std::sqrt(1 - cos_theta * cos_theta);

        Float phi = sample_phi(cos_theta, wi_u, sampler, si);

        assert(phi == phi);
        assert(cos_theta == cos_theta);
        
        vec3 wh = vec3(cos(phi) * sin_theta, sin(phi) * sin_theta, cos_theta);
        vec2 rot = glm::normalize(vec2(wi_u));
        wh = vec3(wh.x * rot.x + wh.y * rot.y, wh.x * rot.y - wh.y * rot.x, wh.z);

        assert(cos_theta > 0.);
        assert(wh.x == wh.x);
        return wh;
    }

    Float MicrograinMicrosurface::lambda(const Float& tau_0, const vec3& wi_u)
    {
        constexpr Float beta = 1.;
        constexpr Float beta2 = 1.;
        constexpr Float beta4 = 1.;
        constexpr Float beta6 = 1.;
        constexpr Float beta8 = 1.;
        constexpr Float beta10 = 1.;
        constexpr Float beta12 = 1.;
        
        Float pi_gamma = -std::log(1. - tau_0);
        Float exp_pi_gamma_minus_one = std::exp(pi_gamma) - 1.;

        Float cos_theta = glm::clamp(wi_u.z, 0.00001f, 0.99999f);
        Float mu = cos_theta / std::sqrt(1. - cos_theta * cos_theta);

        Float mu2 = mu * mu;
        Float mu4 = mu2 * mu2;
        Float mu6 = mu4 * mu2;
        Float mu8 = mu6 * mu2;
        Float mu10 = mu8 * mu2;
        Float mu12 = mu10 * mu2;

        Float beta2_mu2 = beta2 + mu2;
        Float sqrt_beta2_mu2 = std::sqrt(beta2_mu2);

        Float F0 = std::pow(pi_gamma, 1) * (-mu + sqrt_beta2_mu2) / (2 * mu);
        
        Float F1 = std::pow(pi_gamma, 2) * (beta2 + 2 * mu * (mu - sqrt_beta2_mu2)) / (8 * mu * sqrt_beta2_mu2);

        Float F2 = std::pow(pi_gamma, 3) * (3 * beta4 + 12 * beta2 * mu2 + 8 * mu4 - 8 * mu * std::pow(beta2_mu2,3. / 2.)) / (96 * mu * std::pow(beta2_mu2, 3. / 2.));

        Float F3 = std::pow(pi_gamma, 4) * (5 * beta6 + 30 * beta4 * mu2 + 40 * beta2 * mu4 + 16 * mu6 - 16 * mu * std::pow(beta2_mu2, 5./ 2.)) / (768 * mu * std::pow(beta2_mu2, 5. / 2.));

        Float F4 = std::pow(pi_gamma, 5) * (35 * beta8 + 280 * beta6 * mu2 + 560 * beta4 * mu4 + 448 * beta2 * mu6 + 128 * mu8 - 128 * mu * std::pow(beta2_mu2, 7. / 2.)) / (30720 * mu * std::pow(beta2_mu2, 7. / 2.));
            
        Float F5 = std::pow(pi_gamma, 6) * (63 * beta10 + 630 * beta8 * mu2 + 1680 * beta6 * mu4 + 2016 * beta4 * mu6 + 1152 * beta2 * mu8 + 256 * mu10 - 256 * mu * std::pow(beta2_mu2, 9. / 2.)) / (368640 * mu * std::pow(beta2_mu2, 9. / 2.));

        Float F6 = std::pow(pi_gamma, 7) * (231 * beta12 + 2772 * beta10 * mu2 + 9240 * beta8 * mu4 + 14784 * beta6 * mu6 + 12672 * beta4 * mu8 + 5632 * beta2 * mu10 + 1024 * mu12 - 1024 * mu * std::pow(beta2_mu2, 11. / 2.)) / (10321920 * mu * std::pow(beta2_mu2, 11. / 2.));

        return (F0 + F1 + F2 + F3 + F4 + F5 + F6) / exp_pi_gamma_minus_one;
    }

    Float MicrograinMicrosurface::sigma(const vec3& wi_u) 
    {
        return pi * 0.5 * (1. + 1. / glm::clamp(wi_u.z,0.00001f,0.99999f));
    }
    
    Float MicrograinMicrosurface::w_plus(const Float& tau_0, const vec3& wi_u, const vec3& wo_u) {
        Float si_plus_so_minus_sn = pi * 0.5 * (1 / glm::clamp(wi_u.z, 0.00001f, 0.99999f) + 1 / glm::clamp(wo_u.z, 0.00001f, 0.99999f));
        return 1 - std::exp(std::log(1. - tau_0) / pi  * si_plus_so_minus_sn);
    }

    Float MicrograinMicrosurface::sigma_base(const vec3& wh_u) 
    {
        return pi * (1 - wh_u[2] * wh_u[2]);
    }

    Float MicrograinMicrosurface::sigma_shadow(const vec3& wh_u, const vec3& wi_u)
    {
        if (wi_u.z >= 0.99999)
            return 0.;

        Float cos_theta_i_sqr = wi_u.z * wi_u.z;
        Float sin_theta_i_sqr = 1. - cos_theta_i_sqr;
        if (wh_u.z >= std::sqrt(sin_theta_i_sqr))
            return 0.;

        Float tan_theta_i_sqr = sin_theta_i_sqr / cos_theta_i_sqr;

        Float cos_theta_m_sqr = wh_u.z * wh_u.z;
        Float sin_theta_m_sqr = 1. - cos_theta_m_sqr;
        Float tan_theta_m_sqr = sin_theta_m_sqr / cos_theta_m_sqr;

        Float cos_phi_q = -1. / std::sqrt(tan_theta_i_sqr * tan_theta_m_sqr);

        Float q_x = std::sqrt(sin_theta_m_sqr) * cos_phi_q;
        Float q_y = std::sqrt(sin_theta_m_sqr * (1. - cos_phi_q * cos_phi_q));
        Float e_x = std::sqrt(tan_theta_i_sqr) * wh_u.z;
        Float sigma_e = std::acos(glm::clamp(-(q_x - e_x) * wi_u.z,0.f,1.f)) / wi_u.z;
        Float sigma_t = std::abs(e_x * q_y);
        Float sigma_c = (pi - std::acos(cos_phi_q)) * sin_theta_m_sqr;

        return std::max(sigma_e - sigma_t - sigma_c,0.f);
    }

    inline glm::vec2 shadow_intersection_point(const vec3& wh_u, const vec3& wi_u, const vec3& wo_u, const vec3& t, const vec3& b) {
        Float a0 = -t[2] / b[2];
        Float a1 = wh_u[2] / b[2];
        Float a0_sqr = a0 * a0;

        Float d = glm::dot(b, wo_u);
        Float sy = 1. / glm::sqrt(1. - glm::clamp(d*d, 0.00001f, 0.99999f));
        Float sy_sqr = sy * sy;

        Float x = (-a0 * a1 / sy_sqr - glm::sqrt((a0_sqr - a1 * a1) / sy_sqr + 1.)) / (a0_sqr / sy_sqr + 1.);
        Float y = a0 * x + a1;

        vec3 p = t * x + b * y;
        return glm::vec2(p[0], p[1]);
    }

    inline glm::vec2 shadow_intersection_point_0(const vec3& wi, const vec3& wo, const vec3& t, const vec3& b) {
        Float a0 = -t[2] / b[2];
        Float a0_sqr = a0 * a0;

        Float d = glm::dot(b, wo);
        Float sy = 1. / glm::sqrt(1. - glm::clamp(d * d, 0.00001f, 0.99999f));
        Float sy_sqr = sy * sy;

        Float x = (- glm::sqrt(a0_sqr / sy_sqr + 1.)) / (a0_sqr / sy_sqr + 1.);
        Float y = a0 * x;

        vec3 p = t * x + b * y;
        return glm::vec2(p[0], p[1]);
    }

    inline Float area_triangle(const vec2& a, const vec2& b, const vec2& c) {
        Float area = a[0] * (b[1] - c[1]) + b[0] * (c[1] - a[1]) + c[0] * (a[1] - b[1]);
        return 0.5 * std::abs(area);
    }


    inline Float area_sector(const vec2& p1, const vec2& p2, const Float& rad) {
        Float cos_theta = glm::clamp(glm::dot(p1, p2) / (glm::length(p1) * glm::length(p2)), -1.f, 1.f);
        return (std::acos(cos_theta) * rad * rad - std::abs(p1[0] * p2[1] - p2[0] * p1[1])) * 0.5;
    }

    inline Float area_sector_i(const vec2& p, const vec2& pi, const Float& cos_theta_u, const vec3& wi_u) {
        vec2 direction = glm::normalize(vec2(wi_u[0], wi_u[1]));

        Float cos_theta_i = glm::max(wi_u[2], 0.0001f);
        Float tan_theta_i = glm::sqrt(1 - cos_theta_i * cos_theta_i) / cos_theta_i;
        vec2 centered_p = p - direction * tan_theta_i * cos_theta_u;
        vec2 centered_pi = pi - direction * tan_theta_i * cos_theta_u;

        glm::mat2 E = glm::mat2(direction / cos_theta_i, vec2(-direction[1], direction[0]));
        glm::mat2 inv_E = glm::inverse(E);

        centered_p = inv_E * centered_p;
        centered_pi = inv_E * centered_pi;

        return area_sector(centered_p, centered_pi, 1) * glm::determinant(E);
    }

    inline void silhouette_points(const vec3& wh_u, const vec3 wi_u, const vec3& wo_u, vec2& pi_plus, vec2& pi_minus) {
        Float cos_th_i_sqr = wi_u[2]* wi_u[2];
        Float sin_th_i_sqr = 1. - cos_th_i_sqr;
        Float cos_th_m_sqr = wh_u[2] * wh_u[2];
        Float tan_th_i = std::sqrt(sin_th_i_sqr/ cos_th_i_sqr);
        Float sign = glm::sign(wi_u[0] * wo_u[1] - wo_u[0] * wi_u[1]);

        Float x_qi = std::sqrt(1 - cos_th_m_sqr / sin_th_i_sqr) * sign;
        Float y_qi = -wh_u[2] / tan_th_i;
        
        vec2 wi_u_proj = glm::normalize(vec2(wi_u.x, wi_u.y));
        Float sin_ph_i = wi_u_proj.y;
        Float cos_ph_i = wi_u_proj.x;
        pi_plus  = vec2( x_qi * sin_ph_i + y_qi * cos_ph_i, -x_qi * cos_ph_i + y_qi * sin_ph_i);
        pi_minus = vec2(-x_qi * sin_ph_i + y_qi * cos_ph_i,  x_qi * cos_ph_i + y_qi * sin_ph_i);
    }

    inline void silhouette_points_0(const vec3 wi_u, const vec3& wo_u, vec2& pi_plus, vec2& pi_minus) {
        Float x_qi = glm::sign(wi_u[0] * wo_u[1] - wo_u[0] * wi_u[1]);

        vec2 wi_u_proj = glm::normalize(vec2(wi_u.x, wi_u.y));
        Float sin_ph_i = wi_u_proj.y;
        Float cos_ph_i = wi_u_proj.x;
        pi_plus  = vec2(x_qi * sin_ph_i, -x_qi * cos_ph_i);
        pi_minus = vec2(-x_qi * sin_ph_i, x_qi * cos_ph_i);
    }

    inline bool is_left(const vec2& a, const vec2& b, const vec2& c) {
        return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0]) > 0;
    }

    Float MicrograinMicrosurface::sigma_shadow_inter(const vec3& wh_u, const vec3& wi_u, const vec3& wo_u, const vec2& p, const vec2& pi, const vec2& po)
    {
        Float s_t = area_triangle(p, pi, po);
        Float s_i = area_sector_i(p, pi, wh_u[2], wi_u);
        Float s_o = area_sector_i(p, po, wh_u[2], wo_u);
        Float s_c = area_sector(pi, po, std::sqrt(1. - wh_u[2] * wh_u[2]));

        return s_t + s_i + s_o - s_c;
    }

    Float MicrograinMicrosurface::sigma_shadow_inter_0(const vec3& wi_u, const vec3& wo_u, const vec2& p, const vec2& pi, const vec2& po)
    {
        Float s_t = area_triangle(p, pi, po);
        Float s_c = area_sector(pi, po, 1.f);
        Float s_i = area_sector_i(p, pi, 0.f, wi_u);
        Float s_o = area_sector_i(p, po, 0.f, wo_u);

        return s_t + s_i + s_o - s_c;  
    }

    Float MicrograinMicrosurface::sigma_shadow_0(const vec3& wi_u)
    {
        // sigma - sigma_base
        if (wi_u.z >= 0.99999)
            return 0.;
        return pi * 0.5 * (1. / glm::clamp(wi_u.z, 0.00001f, 0.99999f) - 1.);
    }

    Float MicrograinMicrosurface::sigma_inv(const Float& Xi, const vec3& wi_u)
    {
        Float cos_theta_i_sqr = wi_u.z * wi_u.z;
        if (Xi <= pi * cos_theta_i_sqr) {
            return std::sqrt(1. - Xi / pi);
        }

        int idx_theta_i = int(float(res_theta_i) * std::acos(wi_u.z) / (pi * 0.5));
        
        //float pos_Xi = float(res_sigma-1) * Xi / sigma(wi_u);
        //int idx_Xi_1 = std::min(int(pos_Xi) , res_sigma - 1);
        //int idx_Xi_2 = std::min(idx_Xi_1 + 1, res_sigma - 1);
        //float t = pos_Xi - float(idx_Xi_1);
        //float theta_1 = inv_tab_sigma[idx_theta_i][idx_Xi_1];
        //float theta_2 = inv_tab_sigma[idx_theta_i][idx_Xi_2];
        //float theta = (1. - t) * theta_1 + t * theta_2;
        //return  std::cos(theta);

        int idx_theta_m = binary_search<float>( tab_sigma[idx_theta_i], Xi, res_theta_m);
        Float cos_t = std::cos(thetas_m[idx_theta_m]);
        return cos_t;
    }

    Float MicrograinMicrosurface::G1(const vec3& wh_u, const vec3& wi_u, const SurfaceInteraction& si)
    {
        Float tau = tau_0->eval(si);
        if(use_smith || sig_asia_2023)
            return 1. / (1. + lambda(tau,wi_u));
    
        return one_to_many(tau,sigma_shadow(wh_u, wi_u));
        //return std::exp(std::log(1. - tau_0) / pi * sigma_shadow(wh_u, wi_u));
    }

    
    Float MicrograinMicrosurface::G1_0(const vec3& wi_u, const SurfaceInteraction& si)
    {
        Float tau = tau_0->eval(si);
        return one_to_many(tau, sigma_shadow_0(wi_u));
        //return std::exp(std::log(1. - tau_0) / pi * sigma_shadow_0(wi_u));
    }


    Float MicrograinMicrosurface::G2(const vec3& wh_u, const vec3& wi_u, const vec3& wo_u, const SurfaceInteraction& si)
    {
        Float tau = tau_0->eval(si);

        if(use_smith || sig_asia_2023)
            return G1(wh_u, wi_u, si) * G1(wh_u, wo_u, si);
        
        
        Float s_i = sigma_shadow(wh_u, wi_u);
        Float s_o = sigma_shadow(wh_u, wo_u);
        Float s_shadow = s_i + s_o;
        
        if (height_and_direction && s_i > 0. && s_o > 0.) {

            // Compute orthogonal frame
            vec3 b = glm::normalize(wi_u + wo_u);
            vec3 t = glm::normalize(glm::cross(wi_u, wo_u));
            t = t * glm::sign(-t[2]);
            vec3 n = glm::normalize(glm::cross(t, b));

            // Compute silhouette points
            // for i
            vec2 pi_plus;
            vec2 pi_minus;
            silhouette_points(wh_u, wi_u, wo_u, pi_plus, pi_minus);
            // for o
            vec2 po_plus;
            vec2 po_minus;
            silhouette_points(wh_u, wo_u, wi_u, po_plus, po_minus);

            // Compute split line
            vec2 p1 = vec2(-n[2] / n[0] * wh_u[2], 0.);
            vec2 p2 = p1 + vec2(-n[1] / n[0], 1.);

            // Get relative position of all 4 silhouette points
            bool il_pi_p = is_left(p1, p2, pi_plus);
            bool il_pi_m = is_left(p1, p2, pi_minus);
            bool il_po_p = is_left(p1, p2, po_plus);
            bool il_po_m = is_left(p1, p2, po_minus);

            // We should always find : same_side_pi = same_side_po
            bool same_side_pi = il_pi_p == il_pi_m;
            bool same_side_po = il_po_p == il_po_m;
            bool same_side_pio = il_pi_p == il_po_m;
            // Find where a shadow is inside the other
            bool full = same_side_pi && same_side_po && same_side_pio;
            // Find where the shadows are not overlaped
            bool null = same_side_pi && same_side_po && (!same_side_pio);

            if (!null) {
    
                if (full){
                    s_shadow -= std::min(s_i, s_o);
                }
                
                else {
                    // Compute shadow point
                    vec2 p = shadow_intersection_point(wh_u, wi_u, wo_u, t, b);
                    float s_shadow_inter = sigma_shadow_inter(wh_u, wi_u, wo_u, p, pi_plus, po_plus);
                    s_shadow -= s_shadow_inter;
                }

            }            
            

        }
        s_shadow = glm::max(s_shadow, 0.f);
        float g2 = one_to_many(tau, s_shadow);
        assert(g2 == g2);
        assert(g2 <= 1.);
        assert(g2 >= 0.);

        return g2;
        //return std::exp(std::log(1. - tau_0) / pi * s_shadow);
    }

    Float MicrograinMicrosurface::G2_0(const vec3& wi_u, const vec3& wo_u, const SurfaceInteraction& si)
    {
        Float tau = tau_0->eval(si);

        if (use_smith || sig_asia_2023)
            return G1_0(wi_u, si) * G1_0(wo_u, si);


        Float s_i = sigma_shadow_0(wi_u);
        Float s_o = sigma_shadow_0(wo_u);
        Float s_shadow = s_i + s_o;

        if (height_and_direction && s_i > 0. && s_o > 0.) {

            // Compute orthogonal frame
            vec3 b = glm::normalize(wi_u + wo_u);
            vec3 t = glm::normalize(glm::cross(wi_u, wo_u));
            t = t * glm::sign(-t[2]);
            vec3 n = glm::normalize(glm::cross(t, b));

            // Compute silhouette points
            // for i
            vec2 pi_plus;
            vec2 pi_minus;
            silhouette_points_0(wi_u, wo_u, pi_plus, pi_minus);
            // for o
            vec2 po_plus;
            vec2 po_minus;
            silhouette_points_0(wo_u, wi_u, po_plus, po_minus);

           
            Float cos_theta_d = glm::dot(glm::normalize(vec2(wi_u)), glm::normalize(vec2(wo_u)));
            bool colinear_same_side = cos_theta_d > 0.99999f;
            bool colinear_diff_side = cos_theta_d < -0.99999f;

            if (colinear_same_side) {
                s_shadow -= std::min(s_i, s_o);
            }
            else if (!colinear_diff_side) {
                vec2 p = shadow_intersection_point_0(wi_u, wo_u, t, b);
                s_shadow -= sigma_shadow_inter_0(wi_u, wo_u, p, pi_plus, po_plus);
            }

        }
        Float G2_0_ = one_to_many( tau, s_shadow);
        //Float G2_0_ = std::exp(std::log(1. - tau_0) / pi * s_shadow);
        assert(G2_0_ == G2_0_);
        return G2_0_;
    }



}