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

/**
 * @file bssrdf.hpp
 * @brief BSSRDF — サブサーフェス散乱 (dipole 近似)
 *
 * Jensen et al. (2001) "A Practical Model for Subsurface Light Transport"
 * に基づく dipole 拡散近似。
 *
 * 皮膚、蝋、大理石、牛乳など半透明の物質内部での光の散乱をシミュレーション。
 * 光が入射した点から離れた位置で出射する効果 (subsurface scattering) を計算。
 *
 * 主な API:
 *   BSSRDFParams<T>         — 媒質パラメータ (σ_a, σ_s', g, η)
 *   bssrdfDipoleRd(r, params) — 距離 r での拡散反射率
 *   bssrdfPresets::skin()     — 人間の皮膚プリセット
 *   bssrdfPresets::marble()   — 大理石プリセット
 *   BSSRDFProfile<T>         — 散乱プロファイルテーブル (事前計算)
 */

#ifndef CALX_IMAGE_RENDER_BSSRDF_HPP
#define CALX_IMAGE_RENDER_BSSRDF_HPP

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

namespace calx::image {

// ========================================================================
// BSSRDF パラメータ
// ========================================================================

/// サブサーフェス散乱の媒質パラメータ
template<typename T>
struct BSSRDFParams {
    RGB<T> sigma_a;    ///< 吸収係数 [mm^{-1}]
    RGB<T> sigma_s;    ///< 散乱係数 (reduced) [mm^{-1}] — σ_s' = σ_s(1-g)
    T eta = T(1.3);    ///< 相対屈折率 (外部/内部)
    T g = T(0);        ///< Henyey-Greenstein 位相関数の非対称パラメータ

    /// σ_t' = σ_a + σ_s' (減衰係数)
    RGB<T> sigma_t() const {
        return {sigma_a.r + sigma_s.r, sigma_a.g + sigma_s.g, sigma_a.b + sigma_s.b};
    }

    /// 有効アルベド α' = σ_s' / σ_t'
    RGB<T> alpha() const {
        auto st = sigma_t();
        return {(st.r > T(0)) ? sigma_s.r / st.r : T(0),
                (st.g > T(0)) ? sigma_s.g / st.g : T(0),
                (st.b > T(0)) ? sigma_s.b / st.b : T(0)};
    }
};

// ========================================================================
// Fresnel 反射率の拡散積分近似
// ========================================================================

namespace detail_bssrdf {

/// Fdr: 内部拡散反射率の近似 (Egan & Hilgeman 1973 / Jensen 2001 式(5))
template<typename T>
T fresnelDiffuseReflectance(T eta) {
    // 近似式: Fdr ≈ -1.440/η² + 0.710/η + 0.668 + 0.0636η
    if (eta >= T(1)) {
        return T(-1.440) / (eta * eta) + T(0.710) / eta + T(0.668) + T(0.0636) * eta;
    } else {
        // η < 1 (光が密な媒質から疎な媒質に入る場合)
        return T(-0.4399) + T(0.7099) / eta - T(0.3319) / (eta * eta)
             + T(0.0636) / (eta * eta * eta);
    }
}

/// A パラメータ: (1 + Fdr) / (1 - Fdr)
template<typename T>
T dipoleA(T eta) {
    T fdr = fresnelDiffuseReflectance(eta);
    return (T(1) + fdr) / (T(1) - fdr);
}

} // namespace detail_bssrdf

// ========================================================================
// Dipole 拡散近似
// ========================================================================

/// 距離 r [mm] での拡散反射率 Rd(r) (チャネルごと)
/// Jensen 2001 式(4): 実光源 + 仮想光源の dipole
template<typename T>
T bssrdfDipoleRdChannel(T r, T sigma_a, T sigma_s_prime, T eta) {
    T sigma_t_prime = sigma_a + sigma_s_prime;
    if (sigma_t_prime <= T(0)) return T(0);

    T alpha_prime = sigma_s_prime / sigma_t_prime;
    T A = detail_bssrdf::dipoleA(eta);

    // 有効輸送係数
    T sigma_tr = std::sqrt(T(3) * sigma_a * sigma_t_prime);

    // 光源深度
    T one_over_sigma_t = T(1) / sigma_t_prime;

    // 実光源の深度 z_r = 1/σ_t'
    T z_r = one_over_sigma_t;

    // 仮想光源の深度 z_v = z_r + 4A·D = z_r(1 + 4A/3)
    // D = 1/(3σ_t') → z_v = 1/σ_t' + 4A/(3σ_t') = (1 + 4A/3)/σ_t'
    T z_v = one_over_sigma_t * (T(1) + T(4) * A / T(3));

    // 実光源からの距離
    T d_r = std::sqrt(r * r + z_r * z_r);
    // 仮想光源からの距離
    T d_v = std::sqrt(r * r + z_v * z_v);

    T inv_4pi = T(1) / (T(4) * T(3.14159265358979323846));

    // Rd(r) = (α'/4π) [z_r(σ_tr d_r + 1)e^{-σ_tr d_r}/d_r³
    //                  + z_v(σ_tr d_v + 1)e^{-σ_tr d_v}/d_v³]
    T term_r = z_r * (sigma_tr * d_r + T(1)) * std::exp(-sigma_tr * d_r)
               / (d_r * d_r * d_r);
    T term_v = z_v * (sigma_tr * d_v + T(1)) * std::exp(-sigma_tr * d_v)
               / (d_v * d_v * d_v);

    return alpha_prime * inv_4pi * (term_r + term_v);
}

/// RGB 版: 各チャネル独立に計算
template<typename T>
RGB<T> bssrdfDipoleRd(T r, const BSSRDFParams<T>& params) {
    return {
        bssrdfDipoleRdChannel(r, params.sigma_a.r, params.sigma_s.r, params.eta),
        bssrdfDipoleRdChannel(r, params.sigma_a.g, params.sigma_s.g, params.eta),
        bssrdfDipoleRdChannel(r, params.sigma_a.b, params.sigma_s.b, params.eta)
    };
}

/// 総拡散反射率 (全距離積分): Rd_total = ∫_0^∞ 2π·r·Rd(r) dr
/// 数値積分版 (台形則)
template<typename T>
RGB<T> bssrdfTotalReflectance(const BSSRDFParams<T>& params,
                               T maxRadius = T(25), int samples = 200) {
    RGB<T> sum{T(0), T(0), T(0)};
    T dr = maxRadius / samples;
    for (int i = 0; i < samples; ++i) {
        T r = (T(i) + T(0.5)) * dr;
        auto rd = bssrdfDipoleRd(r, params);
        T weight = T(2) * T(3.14159265358979323846) * r * dr;
        sum.r += rd.r * weight;
        sum.g += rd.g * weight;
        sum.b += rd.b * weight;
    }
    return sum;
}

// ========================================================================
// 散乱プロファイルテーブル (事前計算)
// ========================================================================

/// 距離ごとの散乱プロファイルを事前計算して保持
template<typename T>
class BSSRDFProfile {
public:
    BSSRDFProfile() = default;

    /// プロファイルを構築
    /// @param params     媒質パラメータ
    /// @param maxRadius  最大サンプリング半径 [mm]
    /// @param numSamples サンプル数
    BSSRDFProfile(const BSSRDFParams<T>& params,
                  T maxRadius = T(25), int numSamples = 256)
        : maxRadius_(maxRadius), dr_(maxRadius / numSamples) {
        table_.resize(numSamples);
        for (int i = 0; i < numSamples; ++i) {
            T r = (T(i) + T(0.5)) * dr_;
            table_[i] = bssrdfDipoleRd(r, params);
        }
    }

    /// 距離 r でのプロファイル値を線形補間で取得
    RGB<T> sample(T r) const {
        if (table_.empty() || r < T(0)) return {T(0), T(0), T(0)};
        T idx = r / dr_ - T(0.5);
        if (idx < T(0)) idx = T(0);
        int i0 = static_cast<int>(idx);
        if (i0 >= static_cast<int>(table_.size()) - 1)
            return table_.back();
        T frac = idx - T(i0);
        return {
            table_[i0].r * (T(1) - frac) + table_[i0 + 1].r * frac,
            table_[i0].g * (T(1) - frac) + table_[i0 + 1].g * frac,
            table_[i0].b * (T(1) - frac) + table_[i0 + 1].b * frac
        };
    }

    T maxRadius() const { return maxRadius_; }
    int size() const { return static_cast<int>(table_.size()); }

private:
    std::vector<RGB<T>> table_;
    T maxRadius_ = T(25);
    T dr_ = T(0);
};

// ========================================================================
// プリセット
// ========================================================================

namespace bssrdfPresets {

/// 人間の皮膚 (Caucasian skin, Jensen 2001 Table 1)
template<typename T>
BSSRDFParams<T> skin() {
    return {
        {T(0.032), T(0.17), T(0.48)},    // σ_a [mm^-1] (R,G,B)
        {T(0.74), T(0.88), T(1.01)},     // σ_s' [mm^-1]
        T(1.3),                            // η
        T(0)                               // g
    };
}

/// 大理石 (Jensen 2001)
template<typename T>
BSSRDFParams<T> marble() {
    return {
        {T(0.0021), T(0.0041), T(0.0071)},
        {T(2.19), T(2.62), T(3.00)},
        T(1.5),
        T(0)
    };
}

/// 牛乳 (whole milk, Jensen 2001)
template<typename T>
BSSRDFParams<T> milk() {
    return {
        {T(0.0015), T(0.0046), T(0.019)},
        {T(2.55), T(3.21), T(3.77)},
        T(1.3),
        T(0)
    };
}

/// 蝋 (wax)
template<typename T>
BSSRDFParams<T> wax() {
    return {
        {T(0.005), T(0.020), T(0.050)},
        {T(1.0), T(1.2), T(1.4)},
        T(1.4),
        T(0)
    };
}

/// ジェイド (翡翠, jade)
template<typename T>
BSSRDFParams<T> jade() {
    return {
        {T(0.015), T(0.030), T(0.080)},
        {T(0.66), T(0.83), T(0.50)},
        T(1.6),
        T(0)
    };
}

/// カスタムパラメータ: 平均自由行程 (mfp) と色から推定
/// @param mfp     平均自由行程 [mm] — 光が散乱せずに進む平均距離
/// @param albedo  表面アルベド (散乱色)
/// @param eta     屈折率
template<typename T>
BSSRDFParams<T> fromMeanFreePath(T mfp, const RGB<T>& albedo, T eta = T(1.3)) {
    // σ_t' ≈ 1/mfp, σ_s' = α'·σ_t', σ_a = (1-α')·σ_t'
    T sigma_t = T(1) / std::max(mfp, T(0.001));
    return {
        {sigma_t * (T(1) - albedo.r), sigma_t * (T(1) - albedo.g), sigma_t * (T(1) - albedo.b)},
        {sigma_t * albedo.r, sigma_t * albedo.g, sigma_t * albedo.b},
        eta,
        T(0)
    };
}

} // namespace bssrdfPresets

} // namespace calx::image

#endif // CALX_IMAGE_RENDER_BSSRDF_HPP
