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

// scene_renderer.hpp — 科学可視化用シーンレンダラー
//
// 数学的な曲面やデータを美しくレンダリングする高レベル API。
// RayTracer + BVH + ColorMap を統合。
//
// 使用例:
//   SceneRenderer<double> scene(1280, 960);
//   scene.addSurface(
//       [](double x, double y) { return sin(x) * cos(y); },
//       {.xRange = {-3, 3}, .yRange = {-3, 3},
//        .colorMap = ColorMap<double>::viridis()}
//   );
//   scene.render();
//   scene.savePng("output.png");

#ifndef CALX_IMAGE_RENDER_SCENE_RENDERER_HPP
#define CALX_IMAGE_RENDER_SCENE_RENDERER_HPP

#include <image/render/raytracer.hpp>
#include <image/render/bvh.hpp>
#include <image/types/colormap.hpp>
#include <image/io/png.hpp>
#include <image/io/bmp.hpp>

#include <array>
#include <cmath>
#include <functional>
#include <limits>
#include <vector>

namespace calx::image {

// ========================================================================
// 曲面オプション
// ========================================================================

template<typename T>
struct SurfaceOptions {
    std::array<T, 2> xRange{T(-1), T(1)};
    std::array<T, 2> yRange{T(-1), T(1)};
    int gridN = 100;                         // グリッド分割数
    ColorMap<T> colorMap = ColorMap<T>::viridis();
    T opacity = T(1);
    bool autoRange = true;                   // 色の範囲を自動設定
    T vmin = T(0), vmax = T(1);             // 手動色範囲

    // --- 適応的分割オプション ---
    bool adaptive = false;                   // true で適応的分割を有効化
    T minTriArea = T(1e-4);                  // 最小三角形面積 (停止条件)
    T adaptiveError = T(0.1);                // 高さ誤差の閾値
    int maxAdaptiveDepth = 8;                // 最大再帰深度
};

template<typename T>
struct ParametricOptions {
    std::array<T, 2> uRange{T(0), T(6.283185307179586)};
    std::array<T, 2> vRange{T(0), T(6.283185307179586)};
    int gridN = 100;
    ColorMap<T> colorMap = ColorMap<T>::viridis();
    T opacity = T(1);
    bool closedU = false;    // u 方向が閉じている (シリンダ、トーラス)
    bool closedV = false;    // v 方向が閉じている
};

// ========================================================================
// テッセレーション
// ========================================================================

namespace detail_scene {

    /// z = f(x,y) → 三角形メッシュ
    template<typename T, typename Func>
    std::vector<ColoredTriangle<T>> tessellateSurface(
        Func&& f, const SurfaceOptions<T>& opts)
    {
        int N = opts.gridN;
        T x0 = opts.xRange[0], x1 = opts.xRange[1];
        T y0 = opts.yRange[0], y1 = opts.yRange[1];
        T dx = (x1 - x0) / T(N);
        T dy = (y1 - y0) / T(N);

        // 頂点座標と値の計算
        int vN = N + 1;
        std::vector<Point3D<T>> verts(vN * vN);
        std::vector<T> values(vN * vN);

        T zMin = std::numeric_limits<T>::max();
        T zMax = std::numeric_limits<T>::lowest();

        for (int iy = 0; iy < vN; ++iy) {
            for (int ix = 0; ix < vN; ++ix) {
                T x = x0 + ix * dx;
                T y = y0 + iy * dy;
                T z = f(x, y);
                verts[iy * vN + ix] = {x, z, y};  // y-up: z値をY軸に
                values[iy * vN + ix] = z;
                if (z < zMin) zMin = z;
                if (z > zMax) zMax = z;
            }
        }

        // 色範囲
        T cMin = opts.autoRange ? zMin : opts.vmin;
        T cMax = opts.autoRange ? zMax : opts.vmax;
        if (cMax <= cMin) cMax = cMin + T(1);

        // 頂点カラー
        std::vector<RGB<T>> colors(vN * vN);
        for (int i = 0; i < vN * vN; ++i)
            colors[i] = opts.colorMap.map(values[i], cMin, cMax);

        // 頂点法線 (隣接面の面積重み平均)
        std::vector<Point3D<T>> normals(vN * vN, {T(0), T(0), T(0)});

        auto addFaceNormal = [&](int i0, int i1, int i2) {
            auto e1 = verts[i1] - verts[i0];
            auto e2 = verts[i2] - verts[i0];
            auto n = cross(e1, e2);
            normals[i0] = normals[i0] + n;
            normals[i1] = normals[i1] + n;
            normals[i2] = normals[i2] + n;
        };

        for (int iy = 0; iy < N; ++iy) {
            for (int ix = 0; ix < N; ++ix) {
                int i00 = iy * vN + ix;
                int i10 = iy * vN + ix + 1;
                int i01 = (iy + 1) * vN + ix;
                int i11 = (iy + 1) * vN + ix + 1;
                addFaceNormal(i00, i10, i11);
                addFaceNormal(i00, i11, i01);
            }
        }

        for (auto& n : normals) {
            T len = n.norm();
            if (len > T(1e-10)) n = n * (T(1) / len);
        }

        // 三角形生成
        std::vector<ColoredTriangle<T>> tris;
        tris.reserve(N * N * 2);

        for (int iy = 0; iy < N; ++iy) {
            for (int ix = 0; ix < N; ++ix) {
                int i00 = iy * vN + ix;
                int i10 = iy * vN + ix + 1;
                int i01 = (iy + 1) * vN + ix;
                int i11 = (iy + 1) * vN + ix + 1;

                tris.push_back({
                    verts[i00], verts[i10], verts[i11],
                    colors[i00], colors[i10], colors[i11],
                    normals[i00], normals[i10], normals[i11],
                    opts.opacity
                });
                tris.push_back({
                    verts[i00], verts[i11], verts[i01],
                    colors[i00], colors[i11], colors[i01],
                    normals[i00], normals[i11], normals[i01],
                    opts.opacity
                });
            }
        }

        return tris;
    }

    /// 適応的テッセレーション (balanced subdivision)
    /// Phase 1: 各セルの必要分割深度を計算
    /// Phase 2: 隣接セル間の深度差を 1 以下にバランシング
    /// Phase 3: 共有辺の頂点をスナップして T ジャンクション (隙間) を除去
    template<typename T, typename Func>
    std::vector<ColoredTriangle<T>> tessellateAdaptive(
        Func&& f, const SurfaceOptions<T>& opts)
    {
        int N = opts.gridN;
        T x0 = opts.xRange[0], x1 = opts.xRange[1];
        T y0 = opts.yRange[0], y1 = opts.yRange[1];
        T cellW = (x1 - x0) / T(N);
        T cellH = (y1 - y0) / T(N);
        T minArea = opts.minTriArea;
        T errThresh = opts.adaptiveError;
        int maxDepth = opts.maxAdaptiveDepth;

        T zMin = std::numeric_limits<T>::max();
        T zMax = std::numeric_limits<T>::lowest();
        auto evalF = [&](T x, T y) -> T {
            T z = f(x, y);
            if (z < zMin) zMin = z;
            if (z > zMax) zMax = z;
            return z;
        };

        // ---- Phase 1: 各セルの必要分割深度を計算 ----
        std::vector<std::vector<int>> depthMap(N, std::vector<int>(N, 0));

        std::function<int(T, T, T, T, int)> needDepth =
            [&](T ax, T ay, T bx, T by, int d) -> int {
            if (d >= maxDepth) return d;
            T area = (bx - ax) * (by - ay);
            if (area < minArea) return d;
            T mx = (ax + bx) / T(2), my = (ay + by) / T(2);
            T z00 = evalF(ax, ay), z10 = evalF(bx, ay);
            T z01 = evalF(ax, by), z11 = evalF(bx, by);
            T zm = evalF(mx, my);
            if (std::abs(zm - (z00 + z10 + z01 + z11) / T(4)) < errThresh)
                return d;
            return std::max({
                needDepth(ax, ay, mx, my, d + 1),
                needDepth(mx, ay, bx, my, d + 1),
                needDepth(ax, my, mx, by, d + 1),
                needDepth(mx, my, bx, by, d + 1)});
        };

        for (int j = 0; j < N; ++j)
            for (int i = 0; i < N; ++i)
                depthMap[j][i] = needDepth(
                    x0 + i * cellW, y0 + j * cellH,
                    x0 + (i + 1) * cellW, y0 + (j + 1) * cellH, 0);

        // ---- Phase 2: バランシング (隣接差 ≤ 1) ----
        for (bool changed = true; changed;) {
            changed = false;
            for (int j = 0; j < N; ++j)
                for (int i = 0; i < N; ++i) {
                    int d = depthMap[j][i];
                    auto bump = [&](int nj, int ni) {
                        if (nj >= 0 && nj < N && ni >= 0 && ni < N &&
                            depthMap[nj][ni] < d - 1) {
                            depthMap[nj][ni] = d - 1;
                            changed = true;
                        }
                    };
                    bump(j - 1, i); bump(j + 1, i);
                    bump(j, i - 1); bump(j, i + 1);
                }
        }

        // ---- Phase 3: メッシュ生成 (辺スナップ付き) ----
        struct RawTri {
            Point3D<T> v0, v1, v2;
            T z0, z1, z2;
        };
        std::vector<RawTri> rawTris;

        for (int cj = 0; cj < N; ++cj) {
            for (int ci = 0; ci < N; ++ci) {
                int d = depthMap[cj][ci];
                int subdiv = 1 << d;
                T cx0 = x0 + ci * cellW, cy0 = y0 + cj * cellH;
                T sdx = cellW / T(subdiv), sdy = cellH / T(subdiv);
                int vn = subdiv + 1;

                // サブグリッドの頂点を評価
                std::vector<T> zv(vn * vn);
                std::vector<Point3D<T>> verts(vn * vn);
                for (int sj = 0; sj < vn; ++sj)
                    for (int si = 0; si < vn; ++si) {
                        T px = cx0 + si * sdx, py = cy0 + sj * sdy;
                        T z = evalF(px, py);
                        zv[sj * vn + si] = z;
                        verts[sj * vn + si] = {px, z, py};
                    }

                // 隣接セルの深度
                int dL = (ci > 0)     ? depthMap[cj][ci - 1] : d;
                int dR = (ci < N - 1) ? depthMap[cj][ci + 1] : d;
                int dB = (cj > 0)     ? depthMap[cj - 1][ci] : d;
                int dT = (cj < N - 1) ? depthMap[cj + 1][ci] : d;

                // 辺スナップ: 粗い隣接セルとの共有辺で頂点を線形補間
                auto snapEdge = [&](int dNeighbor,
                                    auto getIdx, int edgeLen) {
                    if (dNeighbor >= d) return;  // 同じか細かい → 不要
                    int step = 1 << (d - dNeighbor);
                    for (int k = 0; k < edgeLen; ++k) {
                        if (k % step != 0) {
                            int lo = (k / step) * step;
                            int hi = std::min(lo + step, edgeLen - 1);
                            T t = T(k - lo) / T(hi - lo);
                            int iLo = getIdx(lo), iHi = getIdx(hi);
                            int iK = getIdx(k);
                            T z = zv[iLo] * (T(1) - t) + zv[iHi] * t;
                            verts[iK].y = z;
                        }
                    }
                };

                // 左辺 (si=0), 右辺 (si=subdiv), 下辺 (sj=0), 上辺 (sj=subdiv)
                snapEdge(dL, [&](int k) { return k * vn; }, vn);
                snapEdge(dR, [&](int k) { return k * vn + subdiv; }, vn);
                snapEdge(dB, [&](int k) { return k; }, vn);
                snapEdge(dT, [&](int k) { return subdiv * vn + k; }, vn);

                // 三角形生成
                for (int sj = 0; sj < subdiv; ++sj)
                    for (int si = 0; si < subdiv; ++si) {
                        int i00 = sj * vn + si;
                        int i10 = i00 + 1;
                        int i01 = i00 + vn;
                        int i11 = i01 + 1;
                        rawTris.push_back({verts[i00], verts[i10], verts[i11],
                                           zv[i00], zv[i10], zv[i11]});
                        rawTris.push_back({verts[i00], verts[i11], verts[i01],
                                           zv[i00], zv[i11], zv[i01]});
                    }
            }
        }

        // ---- 色範囲 ----
        T cMin = opts.autoRange ? zMin : opts.vmin;
        T cMax = opts.autoRange ? zMax : opts.vmax;
        if (cMax <= cMin) cMax = cMin + T(1);

        // ---- RawTri → ColoredTriangle ----
        std::vector<ColoredTriangle<T>> tris;
        tris.reserve(rawTris.size());
        for (const auto& rt : rawTris) {
            auto e1 = rt.v1 - rt.v0;
            auto e2 = rt.v2 - rt.v0;
            auto n = cross(e1, e2);
            T len = n.norm();
            if (len > T(1e-10)) n = n * (T(1) / len);
            if (n.y < T(0)) n = n * T(-1);

            tris.push_back({
                rt.v0, rt.v1, rt.v2,
                opts.colorMap.map(rt.z0, cMin, cMax),
                opts.colorMap.map(rt.z1, cMin, cMax),
                opts.colorMap.map(rt.z2, cMin, cMax),
                n, n, n,
                opts.opacity
            });
        }
        return tris;
    }

    /// パラメトリック曲面 → 三角形メッシュ
    /// closedU/closedV: 周期的境界で頂点を共有しシームを除去
    template<typename T, typename Func>
    std::vector<ColoredTriangle<T>> tessellateParametric(
        Func&& f, const ParametricOptions<T>& opts)
    {
        int N = opts.gridN;
        T u0 = opts.uRange[0], u1 = opts.uRange[1];
        T v0 = opts.vRange[0], v1 = opts.vRange[1];
        T du = (u1 - u0) / T(N);
        T dv = (v1 - v0) / T(N);

        // closed 方向は最後の列/行を最初と共有 → N 個だけ生成
        int uVerts = opts.closedU ? N : N + 1;
        int vVerts = opts.closedV ? N : N + 1;

        std::vector<Point3D<T>> verts(vVerts * uVerts);
        std::vector<T> values(vVerts * uVerts);

        T zMin = std::numeric_limits<T>::max();
        T zMax = std::numeric_limits<T>::lowest();

        for (int iv = 0; iv < vVerts; ++iv) {
            for (int iu = 0; iu < uVerts; ++iu) {
                T u = u0 + iu * du;
                T v = v0 + iv * dv;
                auto p = f(u, v);
                int idx = iv * uVerts + iu;
                verts[idx] = p;
                values[idx] = p.y;
                if (p.y < zMin) zMin = p.y;
                if (p.y > zMax) zMax = p.y;
            }
        }

        if (zMax <= zMin) zMax = zMin + T(1);

        std::vector<RGB<T>> colors(vVerts * uVerts);
        for (int i = 0; i < vVerts * uVerts; ++i)
            colors[i] = opts.colorMap.map(values[i], zMin, zMax);

        // インデックスのラップ関数
        auto vidx = [&](int iv, int iu) {
            if (opts.closedU) iu = iu % uVerts;
            if (opts.closedV) iv = iv % vVerts;
            return iv * uVerts + iu;
        };

        // 頂点法線 (面法線の累積)
        std::vector<Point3D<T>> normals(vVerts * uVerts, {T(0), T(0), T(0)});

        auto addFN = [&](int i0, int i1, int i2) {
            auto e1 = verts[i1] - verts[i0];
            auto e2 = verts[i2] - verts[i0];
            auto n = cross(e1, e2);
            normals[i0] = normals[i0] + n;
            normals[i1] = normals[i1] + n;
            normals[i2] = normals[i2] + n;
        };

        for (int iv = 0; iv < N; ++iv) {
            for (int iu = 0; iu < N; ++iu) {
                int i00 = vidx(iv, iu);
                int i10 = vidx(iv, iu + 1);
                int i01 = vidx(iv + 1, iu);
                int i11 = vidx(iv + 1, iu + 1);
                addFN(i00, i10, i11);
                addFN(i00, i11, i01);
            }
        }

        for (auto& n : normals) {
            T len = n.norm();
            if (len > T(1e-10)) n = n * (T(1) / len);
        }

        // 三角形生成
        std::vector<ColoredTriangle<T>> tris;
        tris.reserve(N * N * 2);

        for (int iv = 0; iv < N; ++iv) {
            for (int iu = 0; iu < N; ++iu) {
                int i00 = vidx(iv, iu);
                int i10 = vidx(iv, iu + 1);
                int i01 = vidx(iv + 1, iu);
                int i11 = vidx(iv + 1, iu + 1);

                tris.push_back({
                    verts[i00], verts[i10], verts[i11],
                    colors[i00], colors[i10], colors[i11],
                    normals[i00], normals[i10], normals[i11],
                    opts.opacity
                });
                tris.push_back({
                    verts[i00], verts[i11], verts[i01],
                    colors[i00], colors[i11], colors[i01],
                    normals[i00], normals[i11], normals[i01],
                    opts.opacity
                });
            }
        }

        return tris;
    }

} // namespace detail_scene

// ========================================================================
// SceneRenderer
// ========================================================================

template<typename T>
class SceneRenderer {
public:
    SceneRenderer(std::size_t width, std::size_t height)
        : width_(width), height_(height),
          image_(height, width, 3, T(0)) {}

    // --- カメラ ---
    void camera(const Point3D<T>& position,
                const Point3D<T>& lookAt,
                T fovDeg = T(45)) {
        camPos_ = position;
        camLookAt_ = lookAt;
        fov_ = fovDeg;
    }

    // --- 光源 ---
    void addLight(const Point3D<T>& pos, const RGB<T>& color = {T(1), T(1), T(1)}) {
        lights_.push_back({pos, color});
    }

    // --- 背景 ---
    void setBackground(const RGB<T>& top, const RGB<T>& bottom) {
        bgTop_ = top;
        bgBottom_ = bottom;
    }

    // --- 環境光 ---
    void setAmbient(T a) { ambient_ = a; }

    // --- AA ---
    void setAAThreshold(T t) { aaThreshold_ = t; }
    void setMaxAADepth(int d) { maxAADepth_ = d; }

    // --- 曲面追加: z = f(x, y) ---
    template<typename Func>
    void addSurface(Func&& f, const SurfaceOptions<T>& opts = {}) {
        auto tris = opts.adaptive
            ? detail_scene::tessellateAdaptive<T>(std::forward<Func>(f), opts)
            : detail_scene::tessellateSurface<T>(std::forward<Func>(f), opts);
        appendTriangles(std::move(tris));
    }

    // --- パラメトリック曲面: (u,v) → (x,y,z) ---
    template<typename Func>
    void addParametric(Func&& f, const ParametricOptions<T>& opts = {}) {
        auto tris = detail_scene::tessellateParametric<T>(
            std::forward<Func>(f), opts);
        appendTriangles(std::move(tris));
    }

    // --- 三角形メッシュ直接追加 ---
    void addMesh(const std::vector<Point3D<T>>& vertices,
                 const std::vector<std::array<int, 3>>& faces,
                 const std::vector<T>& values = {},
                 const ColorMap<T>& cmap = ColorMap<T>::viridis()) {
        // 値の範囲
        T vmin = std::numeric_limits<T>::max();
        T vmax = std::numeric_limits<T>::lowest();
        if (!values.empty()) {
            for (T v : values) { vmin = std::min(vmin, v); vmax = std::max(vmax, v); }
        }
        if (vmax <= vmin) { vmax = vmin + T(1); }

        // 頂点カラー
        std::vector<RGB<T>> colors(vertices.size(), {T(0.8), T(0.8), T(0.8)});
        if (!values.empty()) {
            for (std::size_t i = 0; i < vertices.size(); ++i)
                colors[i] = cmap.map(values[i], vmin, vmax);
        }

        // 頂点法線
        std::vector<Point3D<T>> normals(vertices.size(), {T(0), T(0), T(0)});
        for (const auto& face : faces) {
            auto e1 = vertices[face[1]] - vertices[face[0]];
            auto e2 = vertices[face[2]] - vertices[face[0]];
            auto n = cross(e1, e2);
            normals[face[0]] = normals[face[0]] + n;
            normals[face[1]] = normals[face[1]] + n;
            normals[face[2]] = normals[face[2]] + n;
        }
        for (auto& n : normals) {
            T len = n.norm();
            if (len > T(1e-10)) n = n * (T(1) / len);
        }

        std::vector<ColoredTriangle<T>> tris;
        tris.reserve(faces.size());
        for (const auto& face : faces) {
            tris.push_back({
                vertices[face[0]], vertices[face[1]], vertices[face[2]],
                colors[face[0]], colors[face[1]], colors[face[2]],
                normals[face[0]], normals[face[1]], normals[face[2]]
            });
        }
        appendTriangles(std::move(tris));
    }

    // --- レンダリング ---
    void render() {
        // BVH 構築
        bvh_.build(std::move(allTriangles_));

        // カメラ座標系
        auto forward = (camLookAt_ - camPos_).normalized();
        auto right = cross(forward, Point3D<T>{T(0), T(1), T(0)}).normalized();
        auto up = cross(right, forward);

        T aspect = T(width_) / T(height_);
        T halfH = std::tan(fov_ * T(3.14159265358979323846) / T(360));
        T halfW = halfH * aspect;

        // デフォルト光源
        if (lights_.empty()) {
            lights_.push_back({{T(5), T(10), T(-5)}, {T(1), T(1), T(1)}});
            lights_.push_back({{T(-3), T(8), T(-8)}, {T(0.4), T(0.4), T(0.4)}});
        }

        #pragma omp parallel for schedule(dynamic, 4)
        for (int py = 0; py < static_cast<int>(height_); ++py) {
            uint32_t rng = (uint32_t(py) * 73856093u) | 1u;
            for (std::size_t px = 0; px < width_; ++px) {
                rng ^= uint32_t(px) * 19349663u;
                auto color = adaptiveSample(
                    T(px), T(py), T(1), 0,
                    forward, right, up, halfW, halfH, rng);
                color = color.clamped();

                image_.at(py, px, 0) = color.r;
                image_.at(py, px, 1) = color.g;
                image_.at(py, px, 2) = color.b;
            }
        }
    }

    void savePng(const std::string& path) const {
        calx::image::savePng(path, image_);
    }

    void saveBmp(const std::string& path) const {
        calx::image::saveBmp(path, image_);
    }

    const Image<T>& getImage() const { return image_; }
    std::size_t triangleCount() const { return allTriangles_.size(); }

    /// 生の ColoredTriangle を直接追加 (座標軸・ラベル等)
    void addRawTriangles(std::vector<ColoredTriangle<T>> tris) {
        appendTriangles(std::move(tris));
    }

private:
    std::size_t width_, height_;
    Image<T> image_;

    Point3D<T> camPos_{T(3), T(3), T(3)};
    Point3D<T> camLookAt_{T(0), T(0), T(0)};
    T fov_ = T(45);

    RGB<T> bgTop_{T(0.95), T(0.97), T(1.0)};
    RGB<T> bgBottom_{T(1.0), T(1.0), T(1.0)};

    T ambient_ = T(0.1);
    T aaThreshold_ = T(0.03);
    int maxAADepth_ = 2;

    std::vector<Light<T>> lights_;
    std::vector<ColoredTriangle<T>> allTriangles_;
    BVH<T> bvh_;
    bool hasTransparency_ = false;

    void appendTriangles(std::vector<ColoredTriangle<T>> tris) {
        for (const auto& tri : tris) {
            if (tri.opacity < T(1) - T(1e-6)) hasTransparency_ = true;
        }
        allTriangles_.insert(allTriangles_.end(),
            std::make_move_iterator(tris.begin()),
            std::make_move_iterator(tris.end()));
    }

    // --- 背景色 ---
    RGB<T> backgroundColor(const Point3D<T>& dir) const {
        T t = T(0.5) * (dir.y + T(1));
        return (T(1) - t) * bgBottom_ + t * bgTop_;
    }

    // --- 1点のシェーディング (Blinn-Phong) ---
    RGB<T> shade(const Point3D<T>& point, const Point3D<T>& dir,
                 const typename BVH<T>::HitInfo& hit) const {
        auto N = hit.normal;
        if (dot(N, dir) > T(0)) N = T(-1) * N;  // 裏面

        RGB<T> result = hit.color * ambient_;  // ambient

        for (const auto& light : lights_) {
            auto toLight = light.position - point;
            T dist = toLight.norm();
            auto L = toLight * (T(1) / dist);

            // シャドウ: 半透明対応 — 透過率を累積
            T transmittance = shadowTransmittance(point, L, dist);
            if (transmittance < T(1e-3)) continue;

            // Diffuse
            T diff = std::max(T(0), dot(N, L));
            result = result + mulRgb(hit.color, light.color) * (T(0.7) * diff * transmittance);

            // Specular
            auto V = T(-1) * dir;
            auto H = (L + V).normalized();
            T spec = std::pow(std::max(T(0), dot(N, H)), T(64));
            result = result + light.color * (T(0.3) * spec * transmittance);
        }

        return result;
    }

    // --- シャドウ透過率 (半透明面を考慮) ---
    T shadowTransmittance(const Point3D<T>& point, const Point3D<T>& L,
                          T dist) const {
        auto shadowOrigin = point + T(1e-4) * L;
        if (!hasTransparency_) {
            // 不透明のみ: 従来の二値シャドウ
            typename BVH<T>::HitInfo shadow;
            return bvh_.intersect(shadowOrigin, L, T(1e-4), dist, shadow) ? T(0) : T(1);
        }

        // 半透明: 全交差点の不透明度を累積
        thread_local std::vector<typename BVH<T>::HitInfo> shadowHits;
        shadowHits.clear();
        bvh_.intersectAll(shadowOrigin, L, T(1e-4), dist, shadowHits);

        T trans = T(1);
        for (const auto& h : shadowHits) {
            trans *= (T(1) - h.opacity);
            if (trans < T(1e-3)) return T(0);
        }
        return trans;
    }

    // --- 1本のレイをトレース ---
    RGB<T> traceRay(const Point3D<T>& origin, const Point3D<T>& dir) const {
        if (!hasTransparency_) {
            // 高速パス: 全オブジェクト不透明
            typename BVH<T>::HitInfo hit;
            if (!bvh_.intersect(origin, dir, T(1e-4), T(1e10), hit))
                return backgroundColor(dir);
            auto point = origin + hit.t * dir;
            return shade(point, dir, hit);
        }

        // 全交差点を列挙し、front-to-back アルファ合成
        thread_local std::vector<typename BVH<T>::HitInfo> hits;
        hits.clear();
        bvh_.intersectAll(origin, dir, T(1e-4), T(1e10), hits);

        if (hits.empty()) return backgroundColor(dir);

        RGB<T> result{T(0), T(0), T(0)};
        T accAlpha = T(0);

        for (const auto& hit : hits) {
            if (accAlpha >= T(0.999)) break;

            auto point = origin + hit.t * dir;
            RGB<T> shaded = shade(point, dir, hit);
            T alpha = hit.opacity;

            // front-to-back compositing: C_out += (1 - α_acc) * α_i * C_i
            T contrib = (T(1) - accAlpha) * alpha;
            result = result + shaded * contrib;
            accAlpha += contrib;
        }

        // 背景の寄与
        if (accAlpha < T(0.999)) {
            result = result + backgroundColor(dir) * (T(1) - accAlpha);
        }

        return result;
    }

    static RGB<T> mulRgb(const RGB<T>& a, const RGB<T>& b) {
        return {a.r * b.r, a.g * b.g, a.b * b.b};
    }

    // --- ピクセル座標 → レイ ---
    RGB<T> sampleAt(T x, T y,
                    const Point3D<T>& forward, const Point3D<T>& right,
                    const Point3D<T>& up, T halfW, T halfH) const {
        T u = (T(2) * x / T(width_)  - T(1)) * halfW;
        T v = (T(1) - T(2) * y / T(height_)) * halfH;
        auto dir = (forward + u * right + v * up).normalized();
        return traceRay(camPos_, dir);
    }

    // --- 適応的 AA ---
    static T pixelRand(uint32_t& state) {
        state ^= state << 13;
        state ^= state >> 17;
        state ^= state << 5;
        return T(state & 0xFFFFu) / T(0x10000u);
    }

    static T colorDiff(const RGB<T>& a, const RGB<T>& b) {
        return std::max({std::abs(a.r - b.r),
                         std::abs(a.g - b.g),
                         std::abs(a.b - b.b)});
    }

    RGB<T> adaptiveSample(T x, T y, T size, int depth,
                          const Point3D<T>& fwd, const Point3D<T>& rt,
                          const Point3D<T>& up, T hw, T hh,
                          uint32_t& rng) const {
        T jit = size * T(0.1);
        auto s = [&](T sx, T sy) {
            return sampleAt(sx, sy, fwd, rt, up, hw, hh);
        };
        RGB<T> c00 = s(x + jit * pixelRand(rng), y + jit * pixelRand(rng));
        RGB<T> c10 = s(x + size + jit * (pixelRand(rng) - T(0.5)),
                       y + jit * pixelRand(rng));
        RGB<T> c01 = s(x + jit * pixelRand(rng),
                       y + size + jit * (pixelRand(rng) - T(0.5)));
        RGB<T> c11 = s(x + size + jit * (pixelRand(rng) - T(0.5)),
                       y + size + jit * (pixelRand(rng) - T(0.5)));
        T half = size * T(0.5);
        RGB<T> cMid = s(x + half + jit * (pixelRand(rng) - T(0.5)),
                        y + half + jit * (pixelRand(rng) - T(0.5)));

        T maxDiff = std::max({colorDiff(c00, cMid), colorDiff(c10, cMid),
                              colorDiff(c01, cMid), colorDiff(c11, cMid)});

        if (maxDiff > aaThreshold_ && depth < maxAADepth_) {
            auto q00 = adaptiveSample(x,        y,        half, depth+1, fwd, rt, up, hw, hh, rng);
            auto q10 = adaptiveSample(x + half, y,        half, depth+1, fwd, rt, up, hw, hh, rng);
            auto q01 = adaptiveSample(x,        y + half, half, depth+1, fwd, rt, up, hw, hh, rng);
            auto q11 = adaptiveSample(x + half, y + half, half, depth+1, fwd, rt, up, hw, hh, rng);
            return (q00 + q10 + q01 + q11) * T(0.25);
        }
        return (c00 + c10 + c01 + c11 + cMid) * T(0.2);
    }
};

} // namespace calx::image

#endif // CALX_IMAGE_RENDER_SCENE_RENDERER_HPP
