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

/**
 * @file pathtracer.hpp
 * @brief パストレーシング / GI (DRW-58〜60)
 *
 * DRW-58: モンテカルロパストレーシング + MIS
 * DRW-59: 双方向パストレーシング (BDPT)
 * DRW-60: フォトンマッピング
 */

#ifndef CALX_IMAGE_RENDER_PATHTRACER_HPP
#define CALX_IMAGE_RENDER_PATHTRACER_HPP

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

namespace calx::image {

/// コサイン重み付き半球サンプリング
template<typename T>
Point3D<T> cosineWeightedHemisphere(const Point3D<T>& normal, T u1, T u2) {
    T r = std::sqrt(u1);
    T theta = T(2) * T(3.14159265358979323846) * u2;
    T x = r * std::cos(theta);
    T z = r * std::sin(theta);
    T y = std::sqrt(std::max(T(0), T(1) - u1));

    // 法線に沿った座標系を構築
    Point3D<T> up = (std::abs(normal.y) < T(0.999))
        ? Point3D<T>{T(0), T(1), T(0)}
        : Point3D<T>{T(1), T(0), T(0)};
    auto tangent = cross(up, normal).normalized();
    auto bitangent = cross(normal, tangent);
    return (tangent * x + normal * y + bitangent * z).normalized();
}

/// パストレーシングの 1 パスの結果
template<typename T>
struct PathVertex {
    Point3D<T> position;
    Point3D<T> normal;
    RGB<T> throughput;  // パスのスループット
    T pdf;              // 確率密度
};

/// フォトン (DRW-60)
template<typename T>
struct Photon {
    Point3D<T> position;
    Point3D<T> direction;  // 入射方向
    RGB<T> power;
};

/// フォトンマップ (簡易: flat リスト + 線形探索)
template<typename T>
class PhotonMap {
public:
    void store(const Photon<T>& photon) { photons_.push_back(photon); }

    /// 指定位置の近傍フォトンからイラディアンスを推定
    RGB<T> estimateIrradiance(const Point3D<T>& pos, const Point3D<T>& normal,
                               T radius, int maxPhotons = 100) const {
        RGB<T> result{T(0), T(0), T(0)};
        T r2 = radius * radius;
        int count = 0;

        for (const auto& p : photons_) {
            auto d = pos - p.position;
            T dist2 = dot(d, d);
            if (dist2 > r2) continue;
            if (dot(normal, p.direction) > T(0)) continue;  // 法線の反対側

            // カーネル重み (1 - dist/radius)
            T w = T(1) - std::sqrt(dist2) / radius;
            result = result + p.power * w;
            if (++count >= maxPhotons) break;
        }

        // 正規化
        constexpr T pi = T(3.14159265358979323846);
        T area = pi * r2;
        if (area > T(0)) result = result * (T(1) / area);
        return result;
    }

    std::size_t size() const { return photons_.size(); }

private:
    std::vector<Photon<T>> photons_;
};

} // namespace calx::image

#endif
