// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

/**
 * @file pbr.hpp
 * @brief 物理ベースレンダリング (DRW-53〜57, MAT-LAYER)
 *
 * DRW-53: GGX/Beckmann マイクロファセット BRDF, Disney Principled BSDF
 * DRW-54: BSSRDF (サブサーフェス散乱の簡易近似)
 * DRW-55: 異方性反射
 * DRW-56: 法線マッピング (接線空間)
 * DRW-57: 発光マテリアル (エミッシブ)
 * MAT-LAYER: LayeredMaterial (Coat + Sheen レイヤー)
 */

#ifndef CALX_IMAGE_RENDER_PBR_HPP
#define CALX_IMAGE_RENDER_PBR_HPP

#include <image/geometry/primitives.hpp>
#include <image/types/RGB.hpp>
#include <cmath>
#include <algorithm>

namespace calx::image {

// ========================================================================
// PBR マテリアル
// ========================================================================

template<typename T>
struct PBRMaterial {
    RGB<T> albedo{T(0.8), T(0.8), T(0.8)};
    T metallic = T(0);       // 0=誘電体, 1=金属
    T roughness = T(0.5);    // 0=鏡面, 1=完全拡散
    T ior = T(1.5);          // 屈折率
    T subsurface = T(0);     // サブサーフェス散乱の強さ (DRW-54)
    T anisotropy = T(0);     // 異方性 [-1, 1] (DRW-55)
    T anisotropyAngle = T(0); // 異方性の回転角 [rad]
    RGB<T> emission{T(0), T(0), T(0)};  // 発光色 (DRW-57)
    T emissionStrength = T(0);

    // 法線マッピング (DRW-56)
    std::function<Point3D<T>(const Point3D<T>&, const Point3D<T>&)> normalMap;
};

// ========================================================================
// BRDF 関数
// ========================================================================

namespace pbr {

    /// Fresnel (Schlick 近似)
    template<typename T>
    RGB<T> fresnelSchlick(T cosTheta, const RGB<T>& F0) {
        T t = std::pow(T(1) - std::clamp(cosTheta, T(0), T(1)), T(5));
        return F0 * (T(1) - t) + RGB<T>{t, t, t};
    }

    /// GGX (Trowbridge-Reitz) 法線分布関数
    template<typename T>
    T distributionGGX(T NdotH, T roughness) {
        T a = roughness * roughness;
        T a2 = a * a;
        T denom = NdotH * NdotH * (a2 - T(1)) + T(1);
        constexpr T pi = T(3.14159265358979323846);
        return a2 / (pi * denom * denom + T(1e-7));
    }

    /// Beckmann 法線分布関数
    template<typename T>
    T distributionBeckmann(T NdotH, T roughness) {
        T a = roughness * roughness;
        T NdotH2 = NdotH * NdotH;
        T t = (NdotH2 - T(1)) / (a * NdotH2 + T(1e-7));
        constexpr T pi = T(3.14159265358979323846);
        return std::exp(t) / (pi * a * NdotH2 * NdotH2 + T(1e-7));
    }

    /// Smith GGX 幾何関数
    template<typename T>
    T geometrySmithGGX(T NdotV, T NdotL, T roughness) {
        auto ggx = [](T NdotX, T k) {
            return NdotX / (NdotX * (T(1) - k) + k + T(1e-7));
        };
        T k = (roughness + T(1)) * (roughness + T(1)) / T(8);
        return ggx(NdotV, k) * ggx(NdotL, k);
    }

    /// Cook-Torrance スペキュラ BRDF
    template<typename T>
    RGB<T> cookTorranceSpecular(const Point3D<T>& N, const Point3D<T>& V,
                                 const Point3D<T>& L, const PBRMaterial<T>& mat) {
        auto H = (V + L).normalized();
        T NdotH = std::max(dot(N, H), T(0));
        T NdotV = std::max(dot(N, V), T(0));
        T NdotL = std::max(dot(N, L), T(0));
        T HdotV = std::max(dot(H, V), T(0));

        // F0: 基本反射率
        T f0 = ((mat.ior - T(1)) / (mat.ior + T(1)));
        f0 = f0 * f0;
        RGB<T> F0 = mat.albedo * mat.metallic + RGB<T>{f0, f0, f0} * (T(1) - mat.metallic);

        T D = distributionGGX(NdotH, mat.roughness);
        T G = geometrySmithGGX(NdotV, NdotL, mat.roughness);
        RGB<T> F = fresnelSchlick(HdotV, F0);

        T denom = T(4) * NdotV * NdotL + T(1e-4);
        return F * (D * G / denom);
    }

    /// PBR シェーディング (Disney Principled 簡易版)
    template<typename T>
    RGB<T> shade(const Point3D<T>& point, const Point3D<T>& normal,
                  const Point3D<T>& viewDir, const Point3D<T>& lightDir,
                  const RGB<T>& lightColor, const PBRMaterial<T>& mat) {
        auto N = normal.normalized();
        auto V = viewDir.normalized();
        auto L = lightDir.normalized();

        // 法線マッピング (DRW-56)
        if (mat.normalMap) N = mat.normalMap(point, N).normalized();

        T NdotL = std::max(dot(N, L), T(0));

        // 拡散: Lambert (金属は拡散なし)
        RGB<T> kD = mat.albedo * (T(1) - mat.metallic);

        // サブサーフェス散乱 (DRW-54: 簡易 wrap lighting)
        if (mat.subsurface > T(0)) {
            T wrap = (dot(N, L) + mat.subsurface) / (T(1) + mat.subsurface);
            wrap = std::max(wrap, T(0));
            NdotL = NdotL * (T(1) - mat.subsurface) + wrap * mat.subsurface;
        }

        // スペキュラ
        RGB<T> spec = cookTorranceSpecular(N, V, L, mat);

        // 合成
        constexpr T pi = T(3.14159265358979323846);
        RGB<T> diffuse = kD / pi;
        RGB<T> result = (diffuse + spec) * lightColor * NdotL;

        // 発光 (DRW-57)
        result = result + mat.emission * mat.emissionStrength;

        return result;
    }

} // namespace pbr

// ========================================================================
// LayeredMaterial — Coat + Sheen レイヤー (MAT-LAYER)
// ========================================================================

/// マルチレイヤーマテリアル
/// ベース PBR の上にクリアコート層とシーン層を追加する。
/// 紙、布、革、ニス塗り木材、車の塗装などの表現に使用。
template<typename T>
struct LayeredMaterial {
    PBRMaterial<T> base;

    // ── テクスチャ (3D 座標 → 色, PBRMaterial の albedo を上書き) ──
    std::function<RGB<T>(const Point3D<T>&)> texture;

    // ── Coat レイヤー (クリアコート: ニス、ラミネート、車の塗装) ──
    T coatWeight    = T(0);     // [0,1] コート強度
    T coatRoughness = T(0.05);  // コート層の粗さ
    T coatIOR       = T(1.5);   // コート層の屈折率
    RGB<T> coatColor{T(1), T(1), T(1)};  // コートのティント

    // ── Sheen レイヤー (縁光り: 布、紙、フェルト) ──
    T sheenWeight    = T(0);    // [0,1] シーン強度
    RGB<T> sheenColor{T(1), T(1), T(1)};  // シーン色
    T sheenRoughness = T(0.5);  // シーン粗さ [0,1]

    // ── Velvet レイヤー (逆フレネル: ベルベット、スエード) ──
    T velvetWeight    = T(0);   // [0,1] ベルベット効果の強度
    RGB<T> velvetColor{T(1), T(1), T(1)};  // ベルベット色
    T velvetRoughness = T(0.5); // ベルベット粗さ [0,1]

    /// 指定点での albedo を取得 (テクスチャがあればそちらを使用)
    RGB<T> albedoAt(const Point3D<T>& p) const {
        return texture ? texture(p) : base.albedo;
    }
};

namespace pbr {

    // ────────────────────────────────────────────────────────────
    // Charlie NDF — Estevez-Kulla sheen モデル (Frostbite/Disney 方式)
    //   sin(θ)^(2/α - 2) を正規化した分布。布・紙の縁光りに使用。
    // ────────────────────────────────────────────────────────────

    /// Charlie (inverted Gaussian) 法線分布 — sheen 用
    template<typename T>
    T distributionCharlie(T NdotH, T roughness) {
        T alpha = std::max(roughness * roughness, T(1e-6));
        T invAlpha = T(1) / alpha;
        constexpr T pi = T(3.14159265358979323846);
        // D = (2 + 1/α) / (2π) · sin(θ)^(1/α)
        // sin²(θ) = 1 - cos²(θ)
        T sinTheta2 = T(1) - NdotH * NdotH;
        T sinTheta = std::sqrt(std::max(sinTheta2, T(0)));
        T D = (T(2) + invAlpha) / (T(2) * pi) *
              std::pow(std::max(sinTheta, T(1e-12)), invAlpha);
        return D;
    }

    /// Sheen の可視性項 (Neubelt の近似)
    template<typename T>
    T visibilitySheen(T NdotV, T NdotL) {
        return T(1) / (T(4) * (NdotL + NdotV - NdotL * NdotV) + T(1e-7));
    }

    /// Sheen BRDF (Charlie NDF + Neubelt V)
    template<typename T>
    RGB<T> sheenBRDF(const Point3D<T>& N, const Point3D<T>& V,
                     const Point3D<T>& L,
                     T sheenWeight, const RGB<T>& sheenColor, T sheenRoughness) {
        if (sheenWeight <= T(0)) return {T(0), T(0), T(0)};

        auto H = (V + L).normalized();
        T NdotH = std::max(dot(N, H), T(0));
        T NdotV = std::max(dot(N, V), T(1e-5));
        T NdotL = std::max(dot(N, L), T(1e-5));

        T D = distributionCharlie(NdotH, sheenRoughness);
        T Vis = visibilitySheen(NdotV, NdotL);

        return sheenColor * (sheenWeight * D * Vis);
    }

    // ────────────────────────────────────────────────────────────
    // Ashikhmin Velvet BRDF (MAT-CLOTH-1)
    //   ベルベット・スエード等の短い起毛繊維による逆フレネル反射。
    //   通常のマテリアルとは逆に、浅い角度 (grazing) で暗く、
    //   縁方向で明るくなる。Ashikhmin & Premoze (2007) のモデル。
    // ────────────────────────────────────────────────────────────

    /// Ashikhmin velvet NDF — 逆ガウス分布
    /// horizonScattering: 縁方向の散乱強度 [0,1]
    template<typename T>
    T distributionVelvet(T NdotH, T roughness) {
        T alpha = std::max(roughness * roughness, T(1e-6));
        T sinTheta2 = T(1) - NdotH * NdotH;
        // 正規化定数を含む逆ガウス: 縁方向 (sinθ→1) でピーク
        T D = (T(1) + T(4) * std::exp(T(-0.5) / (alpha * sinTheta2 + T(1e-8))))
              / (alpha * (T(1) + T(4) * std::exp(T(-0.5) / alpha)));
        return D;
    }

    /// Velvet BRDF (Ashikhmin モデル)
    /// ベルベットの逆フレネル効果: 正面は暗く、縁が明るく光る。
    template<typename T>
    RGB<T> velvetBRDF(const Point3D<T>& N, const Point3D<T>& V,
                      const Point3D<T>& L,
                      T velvetWeight, const RGB<T>& velvetColor, T roughness) {
        if (velvetWeight <= T(0)) return {T(0), T(0), T(0)};

        T NdotV = std::max(dot(N, V), T(1e-5));
        T NdotL = std::max(dot(N, L), T(1e-5));

        // 逆フレネル: 正面で弱く、grazing で強い
        T invFresnelV = T(1) - std::pow(NdotV, T(0.5));
        T invFresnelL = T(1) - std::pow(NdotL, T(0.5));
        T scatter = invFresnelV * invFresnelL;

        // ハーフベクトルベースの分布
        auto H = (V + L).normalized();
        T NdotH = std::max(dot(N, H), T(0));
        T D = distributionVelvet(NdotH, roughness);

        // 正規化
        constexpr T pi = T(3.14159265358979323846);
        T val = velvetWeight * scatter * D / (T(4) * pi);

        return velvetColor * val;
    }

    // ────────────────────────────────────────────────────────────
    // Irawan-Marschner 簡易版 (MAT-CLOTH-3)
    //   糸の断面形状 (楕円) から異方性 BRDF を合成する簡易モデル。
    //   本来は個々の糸のジオメトリを解析するが、ここでは
    //   経糸/緯糸方向の異方性粗さを切り替えて近似する。
    // ────────────────────────────────────────────────────────────

    /// 異方的 GGX NDF (糸方向の異方性を表現)
    template<typename T>
    T distributionGGXAniso(T NdotH, T HdotT, T HdotB, T alphaT, T alphaB) {
        constexpr T pi = T(3.14159265358979323846);
        T a = HdotT / (alphaT + T(1e-7));
        T b = HdotB / (alphaB + T(1e-7));
        T c = NdotH;
        T d = a * a + b * b + c * c;
        return T(1) / (pi * alphaT * alphaB * d * d + T(1e-12));
    }

    /// 糸の微細構造 BRDF (簡易 Irawan-Marschner)
    /// tangent: 糸の方向ベクトル (経糸 or 緯糸)
    /// alphaAlong: 糸方向の粗さ (通常小さい = 滑らか)
    /// alphaAcross: 糸に直交する方向の粗さ (通常大きい = 粗い)
    template<typename T>
    RGB<T> yarnBRDF(const Point3D<T>& N, const Point3D<T>& V,
                    const Point3D<T>& L, const Point3D<T>& tangent,
                    T alphaAlong, T alphaAcross, const RGB<T>& yarnColor) {
        auto H = (V + L).normalized();
        T NdotH = std::max(dot(N, H), T(0));
        T NdotV = std::max(dot(N, V), T(1e-5));
        T NdotL = std::max(dot(N, L), T(1e-5));

        // 接線空間: T (糸方向), B (糸に直交), N (法線)
        auto T_vec = tangent;
        auto B_vec = cross(N, T_vec);
        T bLen = B_vec.norm();
        if (bLen < T(1e-12)) return yarnColor * (T(1) / T(3.14159265));
        B_vec = B_vec * (T(1) / bLen);

        T HdotT = dot(H, T_vec);
        T HdotB = dot(H, B_vec);

        T D = distributionGGXAniso(NdotH, HdotT, HdotB, alphaAlong, alphaAcross);
        T G = geometrySmithGGX(NdotV, NdotL, (alphaAlong + alphaAcross) * T(0.5));

        constexpr T pi = T(3.14159265358979323846);
        T spec = D * G / (T(4) * NdotV * NdotL + T(1e-4));
        // 拡散 + スペキュラの混合
        return yarnColor * (T(1) / pi + spec * T(0.3));
    }

    /// Coat レイヤーの BRDF (GGX スペキュラ)
    template<typename T>
    RGB<T> coatBRDF(const Point3D<T>& N, const Point3D<T>& V,
                    const Point3D<T>& L,
                    T coatWeight, T coatRoughness, T coatIOR,
                    const RGB<T>& coatColor, T& outCoatFresnel) {
        if (coatWeight <= T(0)) {
            outCoatFresnel = T(0);
            return {T(0), T(0), T(0)};
        }

        auto H = (V + L).normalized();
        T NdotH = std::max(dot(N, H), T(0));
        T NdotV = std::max(dot(N, V), T(1e-5));
        T NdotL = std::max(dot(N, L), T(1e-5));
        T HdotV = std::max(dot(H, V), T(0));

        // Coat の F0 (誘電体コート)
        T f0 = (coatIOR - T(1)) / (coatIOR + T(1));
        f0 = f0 * f0;

        // Schlick Fresnel
        T F = f0 + (T(1) - f0) * std::pow(T(1) - HdotV, T(5));
        outCoatFresnel = F * coatWeight;

        T D = distributionGGX(NdotH, coatRoughness);
        T G = geometrySmithGGX(NdotV, NdotL, coatRoughness);

        T denom = T(4) * NdotV * NdotL + T(1e-4);
        T spec = D * G * F / denom;

        return coatColor * (coatWeight * spec);
    }

    // ────────────────────────────────────────────────────────────
    // LayeredMaterial 用シェーディング
    // ────────────────────────────────────────────────────────────

    /// LayeredMaterial のシェーディング
    ///   coat → sheen → base の順にエネルギーを分配する。
    ///   coat のフレネル反射分だけベース+sheen のエネルギーが減る。
    template<typename T>
    RGB<T> shadeLayered(const Point3D<T>& point, const Point3D<T>& normal,
                        const Point3D<T>& viewDir, const Point3D<T>& lightDir,
                        const RGB<T>& lightColor, const LayeredMaterial<T>& mat) {
        auto N = normal.normalized();
        auto V = viewDir.normalized();
        auto L = lightDir.normalized();

        // 法線マッピング
        if (mat.base.normalMap) N = mat.base.normalMap(point, N).normalized();

        T NdotL = std::max(dot(N, L), T(0));

        // SSS (wrap lighting)
        T effectiveNdotL = NdotL;
        if (mat.base.subsurface > T(0)) {
            T wrap = (dot(N, L) + mat.base.subsurface) / (T(1) + mat.base.subsurface);
            wrap = std::max(wrap, T(0));
            effectiveNdotL = NdotL * (T(1) - mat.base.subsurface) + wrap * mat.base.subsurface;
        }

        // ── Coat レイヤー ──
        T coatFresnel = T(0);
        RGB<T> coatContrib = coatBRDF(N, V, L,
            mat.coatWeight, mat.coatRoughness, mat.coatIOR,
            mat.coatColor, coatFresnel);

        // Coat が反射した分だけベースのエネルギーが減る
        T baseAttenuation = T(1) - coatFresnel;

        // ── Sheen レイヤー ──
        RGB<T> sheenContrib = sheenBRDF(N, V, L,
            mat.sheenWeight, mat.sheenColor, mat.sheenRoughness);

        // ── Velvet レイヤー ──
        RGB<T> velvetContrib = velvetBRDF(N, V, L,
            mat.velvetWeight, mat.velvetColor, mat.velvetRoughness);

        // ── Base レイヤー (PBR) ──
        RGB<T> albedo = mat.albedoAt(point);
        RGB<T> kD = albedo * (T(1) - mat.base.metallic);
        // テクスチャ albedo を反映した一時マテリアルでスペキュラ計算
        PBRMaterial<T> baseMat = mat.base;
        baseMat.albedo = albedo;
        RGB<T> spec = cookTorranceSpecular(N, V, L, baseMat);
        constexpr T pi = T(3.14159265358979323846);
        RGB<T> baseContrib = (kD / pi + spec);

        // ── 合成 ──
        // coat はそのまま、sheen + base は coat の透過分だけ
        RGB<T> result = coatContrib * lightColor * NdotL
                      + (baseContrib + sheenContrib + velvetContrib) * baseAttenuation * lightColor * effectiveNdotL;

        // 発光
        result = result + mat.base.emission * mat.base.emissionStrength;

        return result;
    }

} // namespace pbr

} // namespace calx::image

#endif
