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

/**
 * @file raytracer.hpp
 * @brief Whitted-style レイトレーサー
 *
 * calx の既存プリミティブ (Point3D, Sphere, Plane, Triangle3D) と
 * Image<T> + BMP I/O を活用したレイトレーシングエンジン。
 *
 * 対応機能:
 *   - プリミティブ: 球, 平面, 三角形, AABB (箱), トーラス, 二次曲面,
 *     超二次曲面, ディスク, ハイトフィールド, ブロブ, ラス, プリズム,
 *     メッシュ, CSG (union/intersection/difference)
 *   - マテリアル: Phong シェーディング (ambient + diffuse + specular)
 *   - 反射・屈折 (Snell + Fresnel 近似)
 *   - 影 (シャドウレイ)
 *   - 複数光源 (点光源)
 *   - 再帰的レイトレーシング (最大深度指定)
 *
 * 使用例:
 *   using namespace calx::image;
 *   RayTracer<double> rt(800, 600);
 *   rt.camera({0,2,-5}, {0,0,0}, 60);
 *   rt.addLight({{-5,10,-5}, {1,1,1}});
 *   rt.addSphere({{0,1,0}, 1.0}, {.color={1,0,0}, .specular=0.8});
 *   rt.addPlane({{0,0,0}, {0,1,0}}, {.color={0.5,0.5,0.5}});
 *   rt.render();
 *   saveBmp("output.bmp", rt.image());
 */

#ifndef CALX_IMAGE_RENDER_RAYTRACER_HPP
#define CALX_IMAGE_RENDER_RAYTRACER_HPP

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

#include <image/render/texture.hpp>

#include <cmath>
#include <limits>
#include <vector>
#include <variant>
#include <algorithm>
#include <memory>
#include <array>
#include <functional>

namespace calx::image {

// Light が POV.hpp で既に定義されていなければ定義
#ifndef CALX_IMAGE_RENDER_POV_HPP
template<typename T>
struct Light {
    Point3D<T> position;
    RGB<T> color{T(1), T(1), T(1)};
};
#endif

// ========================================================================
// 拡張光源 (POV-16〜19)
// ========================================================================

template<typename T>
struct SceneLight {
    enum class Type { Point, Spot, Area, Directional };
    Type type = Type::Point;

    // 共通
    Point3D<T> position{};
    RGB<T> color{T(1), T(1), T(1)};
    T intensity = T(1);

    // スポットライト (POV-16)
    Point3D<T> direction{T(0), T(-1), T(0)};  // 照射方向
    T falloffAngle = T(45);    // 外側の円錐角 (度)
    T hotspotAngle = T(30);    // 内側の明るい角 (度)
    T tightness = T(1);        // 減衰の鋭さ

    // エリアライト (POV-17)
    Point3D<T> axis1{T(1), T(0), T(0)};  // 面の辺1
    Point3D<T> axis2{T(0), T(0), T(1)};  // 面の辺2
    int samples1 = 4, samples2 = 4;       // サンプル数 (ソフトシャドウ)

    /// スポットライトの減衰係数 [0,1]
    T spotAttenuation(const Point3D<T>& hitPoint) const {
        auto L = (hitPoint - position).normalized();
        auto D = direction.normalized();
        T cosAngle = dot(L, D);
        T cosOuter = std::cos(falloffAngle * T(3.14159265358979323846) / T(180));
        T cosInner = std::cos(hotspotAngle * T(3.14159265358979323846) / T(180));
        if (cosAngle < cosOuter) return T(0);
        if (cosAngle >= cosInner) return T(1);
        T t = (cosAngle - cosOuter) / (cosInner - cosOuter);
        return std::pow(t, tightness);
    }

    /// エリアライトのサンプル点を生成 (ジッタ付き)
    std::vector<Point3D<T>> areaSamples(uint32_t& rng) const {
        std::vector<Point3D<T>> pts;
        pts.reserve(samples1 * samples2);
        for (int i = 0; i < samples1; ++i) {
            for (int j = 0; j < samples2; ++j) {
                // ジッタ付きサンプリング
                auto jitter = [&]() -> T {
                    rng ^= rng << 13; rng ^= rng >> 17; rng ^= rng << 5;
                    return T(rng & 0xFFFFu) / T(0x10000u);
                };
                T u = (T(i) + jitter()) / T(samples1) - T(0.5);
                T v = (T(j) + jitter()) / T(samples2) - T(0.5);
                pts.push_back(position + u * axis1 + v * axis2);
            }
        }
        return pts;
    }

    /// 平行光源の光線方向 (位置に依存しない)
    Point3D<T> directionTo(const Point3D<T>&) const {
        if (type == Type::Directional)
            return (T(-1) * direction).normalized();
        return {};  // 他の型では使わない
    }
};

// ========================================================================
// Ray
// ========================================================================

template<typename T>
struct Ray {
    Point3D<T> origin;
    Point3D<T> direction;  // 正規化済み

    Point3D<T> at(T t) const {
        return origin + t * direction;
    }
};

// ========================================================================
// Material
// ========================================================================

template<typename T>
struct Material {
    RGB<T> color{T(0.8), T(0.8), T(0.8)};
    T ambient   = T(0.1);
    T diffuse   = T(0.7);
    T specular  = T(0.3);
    T shininess = T(32);
    T reflectivity = T(0);     // 0: 反射なし, 1: 完全鏡面
    T transparency = T(0);     // 0: 不透明, 1: 完全透明
    T ior = T(1.5);            // 屈折率 (ガラス ≈ 1.5)

    // --- テクスチャ (POV-11〜14) ---
    std::function<RGB<T>(const Point3D<T>&)> texture;        // 3D 座標 → 色
    std::function<Point3D<T>(const Point3D<T>&,
                              const Point3D<T>&)> bumpMap;   // (点, 法線) → 摂動法線
    bool metallic = false;     // metallic 反射 (反射色 = 表面色)
    T iridescence = T(0);      // 虹色干渉

    /// 指定点での色を取得
    RGB<T> colorAt(const Point3D<T>& p) const {
        return texture ? texture(p) : color;
    }
    /// バンプ適用後の法線
    Point3D<T> normalAt(const Point3D<T>& p, const Point3D<T>& n) const {
        return bumpMap ? bumpMap(p, n) : n;
    }
};

// ========================================================================
// 交差結果
// ========================================================================

template<typename T>
struct HitRecord {
    T t;                       // レイパラメータ
    Point3D<T> point;          // 交差点
    Point3D<T> normal;         // 法線 (外向き)
    const Material<T>* material;
    bool frontFace;            // レイが外側から当たったか
};

// ========================================================================
// シーンオブジェクト (多態なし — variant ベース)
// ========================================================================

template<typename T>
struct RenderSphere {
    Sphere<T> shape;
    Material<T> material;
};

template<typename T>
struct RenderPlane {
    Plane<T> shape;
    Material<T> material;
};

template<typename T>
struct RenderTriangle {
    Triangle3D<T> shape;
    Material<T> material;
};

template<typename T>
struct RenderBox {
    AABB<T> shape;
    Material<T> material;
};

// --- POV-1: トーラス ---
template<typename T>
struct RenderTorus {
    Torus<T> shape;
    Material<T> material;
};

// --- POV-2: 二次曲面 Ax²+By²+Cz²+Dxy+Exz+Fyz+Gx+Hy+Iz+J=0 ---
template<typename T>
struct Quadric {
    T A{}, B{}, C{}, D{}, E{}, F{}, G{}, H{}, I{}, J{};
};
template<typename T>
struct RenderQuadric {
    Quadric<T> shape;
    Material<T> material;
};

// --- POV-3: 超二次曲面 ---
template<typename T>
struct Superellipsoid {
    T e{T(1)};  // 東西パラメータ
    T n{T(1)};  // 南北パラメータ
};
template<typename T>
struct RenderSuperellipsoid {
    Superellipsoid<T> shape;
    Material<T> material;
};

// --- POV-9: ディスク ---
template<typename T>
struct RenderDisc {
    Disc<T> shape;
    Material<T> material;
};

// --- POV-5: ハイトフィールド ---
template<typename T>
struct HeightFieldData {
    std::vector<std::vector<T>> heights;  // heights[row][col]
    AABB<T> bounds;  // ワールド空間のバウンディングボックス
    std::size_t rows() const { return heights.size(); }
    std::size_t cols() const { return heights.empty() ? 0 : heights[0].size(); }
};
template<typename T>
struct RenderHeightField {
    std::shared_ptr<HeightFieldData<T>> shape;
    Material<T> material;
};

// --- POV-6: ブロブ (Metaball) ---
template<typename T>
struct BlobComponent {
    Point3D<T> center;
    T radius{T(1)};
    T strength{T(1)};  // 正=引力, 負=斥力
};
template<typename T>
struct BlobData {
    std::vector<BlobComponent<T>> components;
    T threshold{T(1)};
};
template<typename T>
struct RenderBlob {
    std::shared_ptr<BlobData<T>> shape;
    Material<T> material;
};

// --- POV-7: ラス (回転体) ---
template<typename T>
struct LatheData {
    std::vector<Point2D<T>> points;  // (radius, y) の制御点
};
template<typename T>
struct RenderLathe {
    std::shared_ptr<LatheData<T>> shape;
    Material<T> material;
};

// --- POV-8: プリズム (Y軸方向の押し出し) ---
template<typename T>
struct PrismData {
    std::vector<Point2D<T>> points;  // XZ 平面上の多角形頂点
    T y0{T(0)}, y1{T(1)};           // 押し出し範囲
};
template<typename T>
struct RenderPrism {
    std::shared_ptr<PrismData<T>> shape;
    Material<T> material;
};

// --- POV-10: メッシュ ---
template<typename T>
struct MeshData {
    std::vector<Point3D<T>> vertices;
    std::vector<std::array<int, 3>> faces;      // 三角形インデックス
    std::vector<Point3D<T>> normals;             // per-vertex 法線 (空なら面法線)
    AABB<T> bounds;                              // バウンディングボックス

    void computeBounds() {
        if (vertices.empty()) return;
        bounds.min = bounds.max = vertices[0];
        for (const auto& v : vertices) {
            bounds.min.x = std::min(bounds.min.x, v.x);
            bounds.min.y = std::min(bounds.min.y, v.y);
            bounds.min.z = std::min(bounds.min.z, v.z);
            bounds.max.x = std::max(bounds.max.x, v.x);
            bounds.max.y = std::max(bounds.max.y, v.y);
            bounds.max.z = std::max(bounds.max.z, v.z);
        }
    }
};
template<typename T>
struct RenderMesh {
    std::shared_ptr<MeshData<T>> shape;
    Material<T> material;
};

// --- POV-4: CSG (前方宣言) ---
template<typename T>
struct CsgNode;

template<typename T>
using SceneObject = std::variant<
    RenderSphere<T>,
    RenderPlane<T>,
    RenderTriangle<T>,
    RenderBox<T>,
    RenderTorus<T>,
    RenderQuadric<T>,
    RenderSuperellipsoid<T>,
    RenderDisc<T>,
    RenderHeightField<T>,
    RenderBlob<T>,
    RenderLathe<T>,
    RenderPrism<T>,
    RenderMesh<T>,
    std::shared_ptr<CsgNode<T>>
>;

// --- POV-4: CSG (定義) ---
template<typename T>
struct CsgNode {
    enum class Op { Union, Intersection, Difference };
    Op op;
    std::vector<SceneObject<T>> children;
    Material<T> material;
};

// ========================================================================
// 交差判定
// ========================================================================

namespace detail {

    /// Ray-Sphere 交差
    template<typename T>
    bool intersectSphere(const Ray<T>& ray, const Sphere<T>& s,
                         T tMin, T tMax, T& tHit, Point3D<T>& normal) {
        auto oc = ray.origin - s.center;
        T a = dot(ray.direction, ray.direction);
        T half_b = dot(oc, ray.direction);
        T c = dot(oc, oc) - s.radius * s.radius;
        T disc = half_b * half_b - a * c;
        if (disc < T(0)) return false;

        T sqrtD = std::sqrt(disc);
        T t = (-half_b - sqrtD) / a;
        if (t < tMin || t > tMax) {
            t = (-half_b + sqrtD) / a;
            if (t < tMin || t > tMax) return false;
        }
        tHit = t;
        auto p = ray.at(t);
        normal = (p - s.center) * (T(1) / s.radius);
        return true;
    }

    /// Ray-Plane 交差
    template<typename T>
    bool intersectPlane(const Ray<T>& ray, const Plane<T>& plane,
                        T tMin, T tMax, T& tHit, Point3D<T>& normal) {
        T denom = dot(plane.normal, ray.direction);
        if (std::abs(denom) < T(1e-8)) return false;

        T t = dot(plane.origin - ray.origin, plane.normal) / denom;
        if (t < tMin || t > tMax) return false;
        tHit = t;
        normal = plane.normal;
        return true;
    }

    /// Ray-Triangle 交差 (Möller–Trumbore)
    template<typename T>
    bool intersectTriangle(const Ray<T>& ray, const Triangle3D<T>& tri,
                           T tMin, T tMax, T& tHit, Point3D<T>& normal) {
        auto e1 = tri.p1 - tri.p0;
        auto e2 = tri.p2 - tri.p0;
        auto h = cross(ray.direction, e2);
        T a = dot(e1, h);
        if (std::abs(a) < T(1e-8)) return false;

        T f = T(1) / a;
        auto s = ray.origin - tri.p0;
        T u = f * dot(s, h);
        if (u < T(0) || u > T(1)) return false;

        auto q = cross(s, e1);
        T v = f * dot(ray.direction, q);
        if (v < T(0) || u + v > T(1)) return false;

        T t = f * dot(e2, q);
        if (t < tMin || t > tMax) return false;
        tHit = t;
        normal = cross(e1, e2).normalized();
        return true;
    }

    /// Ray-AABB 交差 (slab method)
    template<typename T>
    bool intersectAABB(const Ray<T>& ray, const AABB<T>& box,
                       T tMin, T tMax, T& tHit, Point3D<T>& normal) {
        T tNear = tMin, tFar = tMax;
        int hitAxis = -1;
        T hitSign = T(1);

        auto check = [&](T origin, T dir, T bmin, T bmax, int axis) {
            if (std::abs(dir) < T(1e-8)) {
                return origin >= bmin && origin <= bmax;
            }
            T inv = T(1) / dir;
            T t0 = (bmin - origin) * inv;
            T t1 = (bmax - origin) * inv;
            T sign = T(-1);
            if (t0 > t1) { std::swap(t0, t1); sign = T(1); }
            if (t0 > tNear) { tNear = t0; hitAxis = axis; hitSign = sign; }
            if (t1 < tFar)  tFar = t1;
            return tNear <= tFar;
        };

        if (!check(ray.origin.x, ray.direction.x, box.min.x, box.max.x, 0)) return false;
        if (!check(ray.origin.y, ray.direction.y, box.min.y, box.max.y, 1)) return false;
        if (!check(ray.origin.z, ray.direction.z, box.min.z, box.max.z, 2)) return false;

        if (tNear < tMin) tNear = tFar;
        if (tNear < tMin || tNear > tMax) return false;

        tHit = tNear;
        normal = Point3D<T>{T(0), T(0), T(0)};
        if (hitAxis == 0) normal.x = hitSign;
        else if (hitAxis == 1) normal.y = hitSign;
        else normal.z = hitSign;
        return true;
    }

    // ==================================================================
    // 多項式ソルバー (トーラス用)
    // ==================================================================

    /// 3次方程式 x³ + ax² + bx + c = 0 の実根 (最大3根、昇順)
    template<typename T>
    int solveCubic(T a, T b, T c, T roots[3]) {
        constexpr T third = T(1) / T(3);
        constexpr T eps = T(1e-10);
        T a2 = a * a;
        // 降次: t = x - a/3 → t³ + pt + q = 0
        T p = b - a2 * third;
        T q = T(2) * a2 * a / T(27) - a * b * third + c;
        T offset = -a * third;

        // 判別式 D = q²/4 + p³/27
        T q2 = q / T(2);
        T p3 = p * third;
        T D = q2 * q2 + p3 * p3 * p3;

        if (D > eps) {
            // 実根1つ
            T sqrtD = std::sqrt(D);
            T u = std::cbrt(-q2 + sqrtD);
            T v = std::cbrt(-q2 - sqrtD);
            roots[0] = u + v + offset;
            return 1;
        } else if (D < -eps) {
            // 実根3つ (三角関数法)
            T mp3 = -p * third;  // positive
            T r = T(2) * std::sqrt(mp3);
            T cosArg = (T(3) * q / (T(2) * p)) * std::sqrt(T(-3) / p);
            cosArg = std::max(T(-1), std::min(T(1), cosArg));
            T phi = std::acos(cosArg) * third;
            constexpr T twoPi3 = T(2.0943951023931953);
            roots[0] = r * std::cos(phi) + offset;
            roots[1] = r * std::cos(phi - twoPi3) + offset;
            roots[2] = r * std::cos(phi - T(2) * twoPi3) + offset;
            if (roots[0] > roots[1]) std::swap(roots[0], roots[1]);
            if (roots[1] > roots[2]) std::swap(roots[1], roots[2]);
            if (roots[0] > roots[1]) std::swap(roots[0], roots[1]);
            return 3;
        } else {
            // D ≈ 0: 重根
            if (std::abs(q2) < eps) {
                roots[0] = offset;
                return 1;
            }
            T u = std::cbrt(-q2);
            roots[0] = T(2) * u + offset;
            roots[1] = -u + offset;
            if (roots[0] > roots[1]) std::swap(roots[0], roots[1]);
            return 2;
        }
    }

    /// 4次方程式 x⁴ + ax³ + bx² + cx + d = 0 の実根 (最大4根、昇順)
    template<typename T>
    int solveQuartic(T a, T b, T c, T d, T roots[4]) {
        // 降次: u = x - a/4 → u⁴ + pu² + qu + r = 0
        T a2 = a * a;
        T p = b - T(3) * a2 / T(8);
        T q = a * a2 / T(8) - a * b / T(2) + c;
        T r = -T(3) * a2 * a2 / T(256) + a2 * b / T(16) - a * c / T(4) + d;
        T offset = -a / T(4);

        if (std::abs(q) < T(1e-10)) {
            // 双二次式: u⁴ + pu² + r = 0
            T disc = p * p - T(4) * r;
            if (disc < T(-1e-10)) return 0;
            if (disc < T(0)) disc = T(0);
            int n = 0;
            T sqrtDisc = std::sqrt(disc);
            T u2_a = (-p + sqrtDisc) / T(2);
            T u2_b = (-p - sqrtDisc) / T(2);
            if (u2_a >= T(-1e-10)) {
                if (u2_a < T(0)) u2_a = T(0);
                T s = std::sqrt(u2_a);
                roots[n++] = -s + offset;
                roots[n++] = s + offset;
            }
            if (u2_b >= T(-1e-10) && std::abs(u2_b - u2_a) > T(1e-10)) {
                if (u2_b < T(0)) u2_b = T(0);
                T s = std::sqrt(u2_b);
                roots[n++] = -s + offset;
                roots[n++] = s + offset;
            }
            std::sort(roots, roots + n);
            return n;
        }

        // Ferrari: 解決3次方程式 y³ - p/2 · y² - r · y + (4pr - q²)/8 = 0
        T cA = -p / T(2);
        T cB = -r;
        T cC = (T(4) * p * r - q * q) / T(8);
        T cubicRoots[3];
        int nCubic = solveCubic(cA, cB, cC, cubicRoots);
        T y = cubicRoots[nCubic - 1];  // 最大実根

        T sqTerm = T(2) * y - p;
        if (sqTerm < T(0)) sqTerm = T(0);
        T w = std::sqrt(sqTerm);

        int n = 0;
        if (w > T(1e-12)) {
            T q2w = q / (T(2) * w);
            // u² + wu + (y + q/(2w)) = 0
            T disc1 = w * w / T(4) - (y + q2w);
            if (disc1 >= T(-1e-10)) {
                if (disc1 < T(0)) disc1 = T(0);
                T s = std::sqrt(disc1);
                roots[n++] = -w / T(2) - s + offset;
                roots[n++] = -w / T(2) + s + offset;
            }
            // u² - wu + (y - q/(2w)) = 0
            T disc2 = w * w / T(4) - (y - q2w);
            if (disc2 >= T(-1e-10)) {
                if (disc2 < T(0)) disc2 = T(0);
                T s = std::sqrt(disc2);
                roots[n++] = w / T(2) - s + offset;
                roots[n++] = w / T(2) + s + offset;
            }
        }
        std::sort(roots, roots + n);
        return n;
    }

    // ==================================================================
    // POV-1: Ray-Torus 交差
    // ==================================================================
    template<typename T>
    bool intersectTorus(const Ray<T>& ray, const Torus<T>& torus,
                        T tMin, T tMax, T& tHit, Point3D<T>& normal) {
        // ローカル座標に変換 (中心を原点へ)
        auto o = ray.origin - torus.center;
        const auto& d = ray.direction;
        T R = torus.majorRadius;
        T r = torus.minorRadius;
        T R2 = R * R, r2 = r * r;

        // |d|²=1 と仮定
        T od = dot(o, d);
        T oo = dot(o, o);
        T m = oo + R2 - r2;

        // 水平成分 (XZ 平面)
        T k = d.x * d.x + d.z * d.z;       // |d_xz|²
        T h = o.x * d.x + o.z * d.z;       // o_xz · d_xz
        T g = o.x * o.x + o.z * o.z;       // |o_xz|²

        // 4次方程式: t⁴ + c3*t³ + c2*t² + c1*t + c0 = 0
        T c3 = T(4) * od;
        T c2 = T(2) * m + c3 * od / T(2) - T(4) * R2 * k;
        // c2 = 2m + 4*od² - 4*R²*k  ... wait let me recompute properly.
        // (t² + 2*od*t + m)² = 4R²(k*t² + 2h*t + g)
        // t⁴ + 4*od*t³ + (4*od² + 2m)*t² + 4*od*m*t + m² = 4R²k*t² + 8R²h*t + 4R²g
        c2 = T(4) * od * od + T(2) * m - T(4) * R2 * k;
        T c1 = T(4) * od * m - T(8) * R2 * h;
        T c0 = m * m - T(4) * R2 * g;

        T roots[4];
        int nRoots = solveQuartic(c3, c2, c1, c0, roots);

        for (int i = 0; i < nRoots; ++i) {
            if (roots[i] >= tMin && roots[i] <= tMax) {
                tHit = roots[i];
                // 法線計算: P 上の最近点 C = R * normalize(P_xz)
                auto p = o + tHit * d;
                T lenXZ = std::sqrt(p.x * p.x + p.z * p.z);
                if (lenXZ < T(1e-12)) lenXZ = T(1e-12);
                Point3D<T> C{R * p.x / lenXZ, T(0), R * p.z / lenXZ};
                normal = (p - C);
                T nLen = normal.norm();
                if (nLen > T(1e-12)) normal = normal * (T(1) / nLen);
                return true;
            }
        }
        return false;
    }

    // ==================================================================
    // POV-2: Ray-Quadric 交差 (Ax²+By²+Cz²+Dxy+Exz+Fyz+Gx+Hy+Iz+J=0)
    // ==================================================================
    template<typename T>
    bool intersectQuadric(const Ray<T>& ray, const Quadric<T>& q,
                          T tMin, T tMax, T& tHit, Point3D<T>& normal) {
        const auto& o = ray.origin;
        const auto& d = ray.direction;

        // f(o+td) = A(ox+t*dx)² + ... は t の2次式 at²+bt+c=0
        T a_coef = q.A * d.x * d.x + q.B * d.y * d.y + q.C * d.z * d.z
                 + q.D * d.x * d.y + q.E * d.x * d.z + q.F * d.y * d.z;

        T b_coef = T(2) * (q.A * o.x * d.x + q.B * o.y * d.y + q.C * o.z * d.z)
                 + q.D * (o.x * d.y + o.y * d.x)
                 + q.E * (o.x * d.z + o.z * d.x)
                 + q.F * (o.y * d.z + o.z * d.y)
                 + q.G * d.x + q.H * d.y + q.I * d.z;

        T c_coef = q.A * o.x * o.x + q.B * o.y * o.y + q.C * o.z * o.z
                 + q.D * o.x * o.y + q.E * o.x * o.z + q.F * o.y * o.z
                 + q.G * o.x + q.H * o.y + q.I * o.z + q.J;

        if (std::abs(a_coef) < T(1e-12)) {
            // 1次式
            if (std::abs(b_coef) < T(1e-12)) return false;
            T t = -c_coef / b_coef;
            if (t < tMin || t > tMax) return false;
            tHit = t;
        } else {
            T disc = b_coef * b_coef - T(4) * a_coef * c_coef;
            if (disc < T(0)) return false;
            T sqrtD = std::sqrt(disc);
            T t1 = (-b_coef - sqrtD) / (T(2) * a_coef);
            T t2 = (-b_coef + sqrtD) / (T(2) * a_coef);
            if (t1 > t2) std::swap(t1, t2);
            if (t1 >= tMin && t1 <= tMax) tHit = t1;
            else if (t2 >= tMin && t2 <= tMax) tHit = t2;
            else return false;
        }

        // 法線: ∇f = (2Ax+Dy+Ez+G, 2By+Dx+Fz+H, 2Cz+Ex+Fy+I)
        auto p = ray.at(tHit);
        normal = {
            T(2) * q.A * p.x + q.D * p.y + q.E * p.z + q.G,
            T(2) * q.B * p.y + q.D * p.x + q.F * p.z + q.H,
            T(2) * q.C * p.z + q.E * p.x + q.F * p.y + q.I
        };
        T nLen = normal.norm();
        if (nLen > T(1e-12)) normal = normal * (T(1) / nLen);
        return true;
    }

    // ==================================================================
    // POV-3: Ray-Superellipsoid 交差 (レイマーチング + Newton)
    // ==================================================================
    template<typename T>
    T superellipsoidF(T x, T y, T z, T e, T n) {
        // f = (|x|^(2/e) + |z|^(2/e))^(e/n) + |y|^(2/n) - 1
        T ie = T(2) / e, in = T(2) / n, en = e / n;
        T xp = std::pow(std::abs(x) + T(1e-30), ie);
        T zp = std::pow(std::abs(z) + T(1e-30), ie);
        T yp = std::pow(std::abs(y) + T(1e-30), in);
        return std::pow(xp + zp + T(1e-30), en) + yp - T(1);
    }

    template<typename T>
    Point3D<T> superellipsoidNormal(T x, T y, T z, T e, T n) {
        // 数値微分
        constexpr T h = T(1e-5);
        T nx = superellipsoidF(x + h, y, z, e, n) - superellipsoidF(x - h, y, z, e, n);
        T ny = superellipsoidF(x, y + h, z, e, n) - superellipsoidF(x, y - h, z, e, n);
        T nz = superellipsoidF(x, y, z + h, e, n) - superellipsoidF(x, y, z - h, e, n);
        T len = std::sqrt(nx * nx + ny * ny + nz * nz);
        if (len < T(1e-12)) return {T(0), T(1), T(0)};
        return {nx / len, ny / len, nz / len};
    }

    template<typename T>
    bool intersectSuperellipsoid(const Ray<T>& ray, const Superellipsoid<T>& se,
                                  T tMin, T tMax, T& tHit, Point3D<T>& normal) {
        // バウンディングボックス [-1,1]³ との交差で範囲を限定
        AABB<T> box{{T(-1), T(-1), T(-1)}, {T(1), T(1), T(1)}};
        T tEnter, tExit;
        Point3D<T> dummy;
        if (!intersectAABB(ray, box, tMin, tMax, tEnter, dummy)) return false;
        // tExit を求める (AABB の出口)
        {
            T tN = tMin, tF = tMax;
            auto check = [&](T orig, T dir, T bmin, T bmax) {
                if (std::abs(dir) < T(1e-8)) return orig >= bmin && orig <= bmax;
                T inv = T(1) / dir;
                T t0 = (bmin - orig) * inv, t1 = (bmax - orig) * inv;
                if (t0 > t1) std::swap(t0, t1);
                if (t0 > tN) tN = t0;
                if (t1 < tF) tF = t1;
                return tN <= tF;
            };
            check(ray.origin.x, ray.direction.x, T(-1), T(1));
            check(ray.origin.y, ray.direction.y, T(-1), T(1));
            check(ray.origin.z, ray.direction.z, T(-1), T(1));
            tExit = tF;
        }
        if (tEnter < tMin) tEnter = tMin;

        // レイマーチング
        constexpr int maxSteps = 128;
        T dt = (tExit - tEnter) / T(maxSteps);
        T prevT = tEnter;
        auto p0 = ray.at(prevT);
        T prevF = superellipsoidF(p0.x, p0.y, p0.z, se.e, se.n);

        for (int i = 1; i <= maxSteps; ++i) {
            T curT = tEnter + dt * T(i);
            auto p = ray.at(curT);
            T curF = superellipsoidF(p.x, p.y, p.z, se.e, se.n);
            if (prevF * curF < T(0)) {
                // 符号変化 → 二分法で精密化
                T lo = prevT, hi = curT;
                for (int j = 0; j < 20; ++j) {
                    T mid = (lo + hi) * T(0.5);
                    auto pm = ray.at(mid);
                    T fMid = superellipsoidF(pm.x, pm.y, pm.z, se.e, se.n);
                    if (prevF * fMid < T(0)) hi = mid;
                    else { lo = mid; prevF = fMid; }
                }
                T t = (lo + hi) * T(0.5);
                if (t >= tMin && t <= tMax) {
                    tHit = t;
                    auto ph = ray.at(t);
                    normal = superellipsoidNormal(ph.x, ph.y, ph.z, se.e, se.n);
                    return true;
                }
            }
            prevT = curT;
            prevF = curF;
        }
        return false;
    }

    // ==================================================================
    // POV-9: Ray-Disc 交差
    // ==================================================================
    template<typename T>
    bool intersectDisc(const Ray<T>& ray, const Disc<T>& disc,
                       T tMin, T tMax, T& tHit, Point3D<T>& normal) {
        // 平面との交差
        T denom = dot(disc.normal, ray.direction);
        if (std::abs(denom) < T(1e-8)) return false;
        T t = dot(disc.center - ray.origin, disc.normal) / denom;
        if (t < tMin || t > tMax) return false;

        // 半径チェック
        auto p = ray.at(t);
        auto v = p - disc.center;
        T dist2 = dot(v, v);
        T outer2 = disc.outerRadius * disc.outerRadius;
        T inner2 = disc.innerRadius * disc.innerRadius;
        if (dist2 > outer2 || dist2 < inner2) return false;

        tHit = t;
        normal = disc.normal;
        return true;
    }

    // ==================================================================
    // POV-5: Ray-HeightField 交差 (DDA グリッド走査)
    // ==================================================================
    template<typename T>
    bool intersectHeightField(const Ray<T>& ray, const HeightFieldData<T>& hf,
                              T tMin, T tMax, T& tHit, Point3D<T>& normal) {
        // 全体 AABB チェック
        T tEntry;
        Point3D<T> entryNorm;
        if (!intersectAABB(ray, hf.bounds, tMin, tMax, tEntry, entryNorm)) return false;

        std::size_t rows = hf.rows(), cols = hf.cols();
        if (rows < 2 || cols < 2) return false;

        // ワールド → グリッド座標変換
        T dx = (hf.bounds.max.x - hf.bounds.min.x) / T(cols - 1);
        T dz = (hf.bounds.max.z - hf.bounds.min.z) / T(rows - 1);
        T yScale = hf.bounds.max.y - hf.bounds.min.y;

        // グリッドセルごとに三角形テスト (簡易版: 全セル走査)
        bool hit = false;
        T closest = tMax;
        for (std::size_t r = 0; r + 1 < rows; ++r) {
            for (std::size_t c = 0; c + 1 < cols; ++c) {
                T x0 = hf.bounds.min.x + T(c) * dx;
                T x1 = x0 + dx;
                T z0 = hf.bounds.min.z + T(r) * dz;
                T z1 = z0 + dz;
                T y00 = hf.bounds.min.y + hf.heights[r][c] * yScale;
                T y10 = hf.bounds.min.y + hf.heights[r][c + 1] * yScale;
                T y01 = hf.bounds.min.y + hf.heights[r + 1][c] * yScale;
                T y11 = hf.bounds.min.y + hf.heights[r + 1][c + 1] * yScale;

                // 2つの三角形
                Triangle3D<T> t1{{x0, y00, z0}, {x1, y10, z0}, {x0, y01, z1}};
                Triangle3D<T> t2{{x1, y10, z0}, {x1, y11, z1}, {x0, y01, z1}};
                T th; Point3D<T> tn;
                if (intersectTriangle(ray, t1, tMin, closest, th, tn)) {
                    closest = th; tHit = th; normal = tn; hit = true;
                }
                if (intersectTriangle(ray, t2, tMin, closest, th, tn)) {
                    closest = th; tHit = th; normal = tn; hit = true;
                }
            }
        }
        return hit;
    }

    // ==================================================================
    // POV-6: Ray-Blob 交差 (レイマーチング)
    // ==================================================================
    template<typename T>
    T blobField(const Point3D<T>& p, const BlobData<T>& blob) {
        T sum = T(0);
        for (const auto& comp : blob.components) {
            auto d = p - comp.center;
            T dist2 = dot(d, d);
            T r2 = comp.radius * comp.radius;
            if (dist2 < r2) {
                T u = T(1) - dist2 / r2;
                sum += comp.strength * u * u;
            }
        }
        return sum;
    }

    template<typename T>
    bool intersectBlob(const Ray<T>& ray, const BlobData<T>& blob,
                       T tMin, T tMax, T& tHit, Point3D<T>& normal) {
        // バウンディングボックスを計算
        if (blob.components.empty()) return false;
        T xMin = T(1e30), yMin = T(1e30), zMin = T(1e30);
        T xMax = T(-1e30), yMax = T(-1e30), zMax = T(-1e30);
        for (const auto& c : blob.components) {
            xMin = std::min(xMin, c.center.x - c.radius);
            yMin = std::min(yMin, c.center.y - c.radius);
            zMin = std::min(zMin, c.center.z - c.radius);
            xMax = std::max(xMax, c.center.x + c.radius);
            yMax = std::max(yMax, c.center.y + c.radius);
            zMax = std::max(zMax, c.center.z + c.radius);
        }
        AABB<T> bb{{xMin, yMin, zMin}, {xMax, yMax, zMax}};
        T tEntry;
        Point3D<T> entryNorm;
        if (!intersectAABB(ray, bb, tMin, tMax, tEntry, entryNorm)) return false;

        // レイマーチング
        T tExit = tMax;
        {
            T tN = tMin, tF = tMax;
            auto check = [&](T orig, T dir, T bmin, T bmax) {
                if (std::abs(dir) < T(1e-8)) return orig >= bmin && orig <= bmax;
                T inv = T(1) / dir;
                T t0 = (bmin - orig) * inv, t1 = (bmax - orig) * inv;
                if (t0 > t1) std::swap(t0, t1);
                if (t0 > tN) tN = t0;
                if (t1 < tF) tF = t1;
                return tN <= tF;
            };
            check(ray.origin.x, ray.direction.x, xMin, xMax);
            check(ray.origin.y, ray.direction.y, yMin, yMax);
            check(ray.origin.z, ray.direction.z, zMin, zMax);
            tExit = tF;
        }

        constexpr int maxSteps = 128;
        T dt = (tExit - tEntry) / T(maxSteps);
        T prevT = tEntry;
        T prevF = blobField(ray.at(prevT), blob) - blob.threshold;

        for (int i = 1; i <= maxSteps; ++i) {
            T curT = tEntry + dt * T(i);
            T curF = blobField(ray.at(curT), blob) - blob.threshold;
            if (prevF * curF < T(0)) {
                // 二分法で精密化
                T lo = prevT, hi = curT;
                T loF = prevF;
                for (int j = 0; j < 20; ++j) {
                    T mid = (lo + hi) * T(0.5);
                    T fMid = blobField(ray.at(mid), blob) - blob.threshold;
                    if (loF * fMid < T(0)) hi = mid;
                    else { lo = mid; loF = fMid; }
                }
                T t = (lo + hi) * T(0.5);
                if (t >= tMin && t <= tMax) {
                    tHit = t;
                    // 数値微分で法線
                    auto p = ray.at(t);
                    constexpr T h = T(1e-4);
                    T fx = blobField({p.x + h, p.y, p.z}, blob)
                         - blobField({p.x - h, p.y, p.z}, blob);
                    T fy = blobField({p.x, p.y + h, p.z}, blob)
                         - blobField({p.x, p.y - h, p.z}, blob);
                    T fz = blobField({p.x, p.y, p.z + h}, blob)
                         - blobField({p.x, p.y, p.z - h}, blob);
                    T nLen = std::sqrt(fx * fx + fy * fy + fz * fz);
                    if (nLen > T(1e-12))
                        normal = {fx / nLen, fy / nLen, fz / nLen};
                    else
                        normal = {T(0), T(1), T(0)};
                    return true;
                }
            }
            prevT = curT;
            prevF = curF;
        }
        return false;
    }

    // ==================================================================
    // POV-7: Ray-Lathe 交差 (Y軸回りの回転体、線形セグメント)
    // ==================================================================
    template<typename T>
    bool intersectLathe(const Ray<T>& ray, const LatheData<T>& lathe,
                        T tMin, T tMax, T& tHit, Point3D<T>& normal) {
        if (lathe.points.size() < 2) return false;
        bool hit = false;
        T closest = tMax;

        // 各線形セグメントの回転による截頭円錐/円柱との交差
        for (std::size_t i = 0; i + 1 < lathe.points.size(); ++i) {
            T r0 = lathe.points[i].x;      // 下端の半径
            T y0 = lathe.points[i].y;      // 下端の Y
            T r1 = lathe.points[i + 1].x;  // 上端の半径
            T y1 = lathe.points[i + 1].y;  // 上端の Y

            // 截頭円錐 f(x,y,z) = x²+z² - R(y)² = 0
            // R(y) = r0 + (r1-r0)*(y-y0)/(y1-y0)
            T dy = y1 - y0;
            if (std::abs(dy) < T(1e-12)) continue;

            T dr = r1 - r0;
            // R(y) = r0 + dr/dy * (y - y0)
            // = A + B*y  where A = r0 - dr*y0/dy, B = dr/dy
            T B = dr / dy;
            T A = r0 - B * y0;

            // x²+z² = (A+B*y)²
            // (ox+t*dx)²+(oz+t*dz)² = (A+B*(oy+t*dy_ray))²
            const auto& o = ray.origin;
            const auto& d = ray.direction;

            T AB = A + B * o.y;
            T Bd = B * d.y;
            // LHS: (dx²+dz²)t² + 2(ox*dx+oz*dz)t + ox²+oz²
            T lhs_a = d.x * d.x + d.z * d.z;
            T lhs_b = T(2) * (o.x * d.x + o.z * d.z);
            T lhs_c = o.x * o.x + o.z * o.z;
            // RHS: (AB+Bd*t)² = Bd²*t² + 2*AB*Bd*t + AB²
            T rhs_a = Bd * Bd;
            T rhs_b = T(2) * AB * Bd;
            T rhs_c = AB * AB;

            T qa = lhs_a - rhs_a;
            T qb = lhs_b - rhs_b;
            T qc = lhs_c - rhs_c;

            T roots[2];
            int nRoots = 0;
            if (std::abs(qa) < T(1e-12)) {
                if (std::abs(qb) > T(1e-12)) {
                    roots[0] = -qc / qb;
                    nRoots = 1;
                }
            } else {
                T disc = qb * qb - T(4) * qa * qc;
                if (disc >= T(0)) {
                    T sqrtD = std::sqrt(disc);
                    roots[0] = (-qb - sqrtD) / (T(2) * qa);
                    roots[1] = (-qb + sqrtD) / (T(2) * qa);
                    nRoots = 2;
                }
            }

            for (int j = 0; j < nRoots; ++j) {
                T t = roots[j];
                if (t < tMin || t > closest) continue;
                T py = o.y + t * d.y;
                T yLo = std::min(y0, y1), yHi = std::max(y0, y1);
                if (py < yLo || py > yHi) continue;
                // 法線
                auto p = ray.at(t);
                T rr = A + B * py;  // R(y)
                T lenXZ = std::sqrt(p.x * p.x + p.z * p.z);
                if (lenXZ < T(1e-12)) continue;
                // 円錐表面の法線: 外向き radial - 傾き分
                Point3D<T> radial{p.x / lenXZ, T(0), p.z / lenXZ};
                Point3D<T> n = radial - B * Point3D<T>{T(0), T(1), T(0)} * (rr > T(0) ? T(1) : T(-1));
                T nLen = n.norm();
                if (nLen > T(1e-12)) n = n * (T(1) / nLen);
                tHit = t; normal = n; closest = t; hit = true;
            }
        }

        // 上端・下端のキャップ (ディスク)
        if (!lathe.points.empty()) {
            T rFirst = lathe.points.front().x;
            T yFirst = lathe.points.front().y;
            T rLast = lathe.points.back().x;
            T yLast = lathe.points.back().y;
            auto testCap = [&](T rCap, T yCap, const Point3D<T>& capNorm) {
                if (std::abs(ray.direction.y) < T(1e-8)) return;
                T t = (yCap - ray.origin.y) / ray.direction.y;
                if (t < tMin || t > closest) return;
                auto p = ray.at(t);
                T dist2 = p.x * p.x + p.z * p.z;
                if (dist2 > rCap * rCap) return;
                tHit = t; normal = capNorm; closest = t; hit = true;
            };
            if (rFirst > T(1e-8))
                testCap(rFirst, yFirst, {T(0), yFirst < yLast ? T(-1) : T(1), T(0)});
            if (rLast > T(1e-8))
                testCap(rLast, yLast, {T(0), yLast > yFirst ? T(1) : T(-1), T(0)});
        }
        return hit;
    }

    // ==================================================================
    // POV-8: Ray-Prism 交差 (Y軸方向の押し出し多角形)
    // ==================================================================
    template<typename T>
    bool intersectPrism(const Ray<T>& ray, const PrismData<T>& prism,
                        T tMin, T tMax, T& tHit, Point3D<T>& normal) {
        if (prism.points.size() < 3) return false;
        bool hit = false;
        T closest = tMax;
        const auto& pts = prism.points;
        std::size_t n = pts.size();

        // 2D 点内包判定 (XZ 平面)
        auto pointInPolygon = [&](T px, T pz) -> bool {
            bool inside = false;
            for (std::size_t i = 0, j = n - 1; i < n; j = i++) {
                T xi = pts[i].x, zi = pts[i].y;  // Point2D.y → Z座標
                T xj = pts[j].x, zj = pts[j].y;
                if ((zi > pz) != (zj > pz) &&
                    px < (xj - xi) * (pz - zi) / (zj - zi) + xi)
                    inside = !inside;
            }
            return inside;
        };

        // 上面・下面キャップ
        auto testCap = [&](T yCap, const Point3D<T>& capNorm) {
            if (std::abs(ray.direction.y) < T(1e-8)) return;
            T t = (yCap - ray.origin.y) / ray.direction.y;
            if (t < tMin || t > closest) return;
            auto p = ray.at(t);
            if (!pointInPolygon(p.x, p.z)) return;
            tHit = t; normal = capNorm; closest = t; hit = true;
        };
        testCap(prism.y0, {T(0), T(-1), T(0)});
        testCap(prism.y1, {T(0), T(1), T(0)});

        // 側面: 各辺の押し出し壁 (無限平面との交差 + 範囲チェック)
        for (std::size_t i = 0; i < n; ++i) {
            std::size_t j = (i + 1) % n;
            T x0 = pts[i].x, z0 = pts[i].y;
            T x1 = pts[j].x, z1 = pts[j].y;
            // 辺の法線 (XZ 平面): (z1-z0, -(x1-x0)) 外向き
            Point3D<T> wallNorm{z1 - z0, T(0), -(x1 - x0)};
            T wLen = std::sqrt(wallNorm.x * wallNorm.x + wallNorm.z * wallNorm.z);
            if (wLen < T(1e-12)) continue;
            wallNorm = wallNorm * (T(1) / wLen);

            // 平面 dot(wallNorm, p - p0) = 0 where p0 = (x0, 0, z0)
            T denom = wallNorm.x * ray.direction.x + wallNorm.z * ray.direction.z;
            if (std::abs(denom) < T(1e-8)) continue;
            T t = (wallNorm.x * (x0 - ray.origin.x) + wallNorm.z * (z0 - ray.origin.z)) / denom;
            if (t < tMin || t > closest) continue;

            auto p = ray.at(t);
            if (p.y < prism.y0 || p.y > prism.y1) continue;

            // 辺上にあるかチェック (パラメータ s ∈ [0,1])
            T ex = x1 - x0, ez = z1 - z0;
            T s;
            if (std::abs(ex) > std::abs(ez))
                s = (p.x - x0) / ex;
            else
                s = (p.z - z0) / ez;
            if (s < T(0) || s > T(1)) continue;

            tHit = t; normal = wallNorm; closest = t; hit = true;
        }
        return hit;
    }

    // ==================================================================
    // POV-10: Ray-Mesh 交差 (全三角形走査)
    // ==================================================================
    template<typename T>
    bool intersectMesh(const Ray<T>& ray, const MeshData<T>& mesh,
                       T tMin, T tMax, T& tHit, Point3D<T>& normal) {
        // AABB による早期棄却
        T bbT; Point3D<T> bbN;
        if (!intersectAABB(ray, mesh.bounds, tMin, tMax, bbT, bbN)) return false;

        bool hit = false;
        T closest = tMax;
        bool hasNormals = (mesh.normals.size() == mesh.vertices.size());

        for (const auto& face : mesh.faces) {
            const auto& v0 = mesh.vertices[face[0]];
            const auto& v1 = mesh.vertices[face[1]];
            const auto& v2 = mesh.vertices[face[2]];

            // Möller–Trumbore
            auto e1 = v1 - v0;
            auto e2 = v2 - v0;
            auto h = cross(ray.direction, e2);
            T a = dot(e1, h);
            if (std::abs(a) < T(1e-8)) continue;
            T f = T(1) / a;
            auto s = ray.origin - v0;
            T u = f * dot(s, h);
            if (u < T(0) || u > T(1)) continue;
            auto q = cross(s, e1);
            T v = f * dot(ray.direction, q);
            if (v < T(0) || u + v > T(1)) continue;
            T t = f * dot(e2, q);
            if (t < tMin || t > closest) continue;

            closest = t;
            tHit = t;
            hit = true;

            if (hasNormals) {
                // 頂点法線の重心補間
                T w = T(1) - u - v;
                normal = (w * mesh.normals[face[0]]
                        + u * mesh.normals[face[1]]
                        + v * mesh.normals[face[2]]).normalized();
            } else {
                normal = cross(e1, e2).normalized();
            }
        }
        return hit;
    }

    // ==================================================================
    // POV-4: CSG 交差 (inside テスト方式)
    // ==================================================================

    /// プリミティブの内側判定
    template<typename T>
    bool isInsideObject(const Point3D<T>& p, const SceneObject<T>& obj) {
        return std::visit([&](const auto& o) -> bool {
            using ObjT = std::decay_t<decltype(o)>;
            if constexpr (std::is_same_v<ObjT, RenderSphere<T>>) {
                auto d = p - o.shape.center;
                return dot(d, d) < o.shape.radius * o.shape.radius;
            } else if constexpr (std::is_same_v<ObjT, RenderPlane<T>>) {
                return dot(p - o.shape.origin, o.shape.normal) < T(0);
            } else if constexpr (std::is_same_v<ObjT, RenderBox<T>>) {
                return o.shape.contains(p);
            } else if constexpr (std::is_same_v<ObjT, RenderTorus<T>>) {
                auto lp = p - o.shape.center;
                T R = o.shape.majorRadius, r = o.shape.minorRadius;
                T S = dot(lp, lp) + R * R - r * r;
                T xz2 = lp.x * lp.x + lp.z * lp.z;
                return S * S - T(4) * R * R * xz2 < T(0);
            } else if constexpr (std::is_same_v<ObjT, RenderQuadric<T>>) {
                const auto& q = o.shape;
                T val = q.A*p.x*p.x + q.B*p.y*p.y + q.C*p.z*p.z
                      + q.D*p.x*p.y + q.E*p.x*p.z + q.F*p.y*p.z
                      + q.G*p.x + q.H*p.y + q.I*p.z + q.J;
                return val < T(0);
            } else if constexpr (std::is_same_v<ObjT, RenderSuperellipsoid<T>>) {
                return superellipsoidF(p.x, p.y, p.z, o.shape.e, o.shape.n) < T(0);
            } else if constexpr (std::is_same_v<ObjT, std::shared_ptr<CsgNode<T>>>) {
                // CSG の内側判定は子要素に委譲
                if (!o) return false;
                const auto& node = *o;
                if (node.op == CsgNode<T>::Op::Union) {
                    for (const auto& child : node.children)
                        if (isInsideObject(p, child)) return true;
                    return false;
                } else if (node.op == CsgNode<T>::Op::Intersection) {
                    for (const auto& child : node.children)
                        if (!isInsideObject(p, child)) return false;
                    return !node.children.empty();
                } else { // Difference
                    if (node.children.empty()) return false;
                    if (!isInsideObject(p, node.children[0])) return false;
                    for (std::size_t i = 1; i < node.children.size(); ++i)
                        if (isInsideObject(p, node.children[i])) return false;
                    return true;
                }
            } else {
                // Triangle, Disc, HeightField, Blob, Lathe, Prism, Mesh
                // — 薄い表面なので内側は常に false
                return false;
            }
        }, obj);
    }

    /// CSG 交差判定: 全子要素の最近交差を調べ、CSG 条件を満たすまで反復
    template<typename T>
    bool intersectCsg(const Ray<T>& ray, const CsgNode<T>& node,
                      T tMin, T tMax, T& tHit, Point3D<T>& normal) {
        if (node.children.empty()) return false;

        if (node.op == CsgNode<T>::Op::Union) {
            // Union: 全子要素の最近交差
            bool hit = false;
            T closest = tMax;
            for (const auto& child : node.children) {
                T th; Point3D<T> tn;
                bool didHit = std::visit([&](const auto& o) -> bool {
                    using ObjT = std::decay_t<decltype(o)>;
                    if constexpr (std::is_same_v<ObjT, RenderSphere<T>>)
                        return intersectSphere(ray, o.shape, tMin, closest, th, tn);
                    else if constexpr (std::is_same_v<ObjT, RenderPlane<T>>)
                        return intersectPlane(ray, o.shape, tMin, closest, th, tn);
                    else if constexpr (std::is_same_v<ObjT, RenderTriangle<T>>)
                        return intersectTriangle(ray, o.shape, tMin, closest, th, tn);
                    else if constexpr (std::is_same_v<ObjT, RenderBox<T>>)
                        return intersectAABB(ray, o.shape, tMin, closest, th, tn);
                    else if constexpr (std::is_same_v<ObjT, RenderTorus<T>>)
                        return intersectTorus(ray, o.shape, tMin, closest, th, tn);
                    else if constexpr (std::is_same_v<ObjT, RenderQuadric<T>>)
                        return intersectQuadric(ray, o.shape, tMin, closest, th, tn);
                    else if constexpr (std::is_same_v<ObjT, RenderSuperellipsoid<T>>)
                        return intersectSuperellipsoid(ray, o.shape, tMin, closest, th, tn);
                    else if constexpr (std::is_same_v<ObjT, RenderDisc<T>>)
                        return intersectDisc(ray, o.shape, tMin, closest, th, tn);
                    else if constexpr (std::is_same_v<ObjT, RenderHeightField<T>>)
                        return o.shape ? intersectHeightField(ray, *o.shape, tMin, closest, th, tn) : false;
                    else if constexpr (std::is_same_v<ObjT, RenderBlob<T>>)
                        return o.shape ? intersectBlob(ray, *o.shape, tMin, closest, th, tn) : false;
                    else if constexpr (std::is_same_v<ObjT, RenderLathe<T>>)
                        return o.shape ? intersectLathe(ray, *o.shape, tMin, closest, th, tn) : false;
                    else if constexpr (std::is_same_v<ObjT, RenderPrism<T>>)
                        return o.shape ? intersectPrism(ray, *o.shape, tMin, closest, th, tn) : false;
                    else if constexpr (std::is_same_v<ObjT, RenderMesh<T>>)
                        return o.shape ? intersectMesh(ray, *o.shape, tMin, closest, th, tn) : false;
                    else if constexpr (std::is_same_v<ObjT, std::shared_ptr<CsgNode<T>>>)
                        return o ? intersectCsg(ray, *o, tMin, closest, th, tn) : false;
                    else
                        return false;
                }, child);
                if (didHit) {
                    closest = th; tHit = th; normal = tn; hit = true;
                }
            }
            return hit;
        }

        // Intersection / Difference: 反復テスト
        constexpr int maxIter = 64;
        T curMin = tMin;
        for (int iter = 0; iter < maxIter; ++iter) {
            // 全子要素の最近交差を探す
            T bestT = tMax + T(1);
            Point3D<T> bestN;
            int bestChild = -1;

            for (int ci = 0; ci < static_cast<int>(node.children.size()); ++ci) {
                T th; Point3D<T> tn;
                bool didHit = std::visit([&](const auto& o) -> bool {
                    using ObjT = std::decay_t<decltype(o)>;
                    if constexpr (std::is_same_v<ObjT, RenderSphere<T>>)
                        return intersectSphere(ray, o.shape, curMin, tMax, th, tn);
                    else if constexpr (std::is_same_v<ObjT, RenderPlane<T>>)
                        return intersectPlane(ray, o.shape, curMin, tMax, th, tn);
                    else if constexpr (std::is_same_v<ObjT, RenderTriangle<T>>)
                        return intersectTriangle(ray, o.shape, curMin, tMax, th, tn);
                    else if constexpr (std::is_same_v<ObjT, RenderBox<T>>)
                        return intersectAABB(ray, o.shape, curMin, tMax, th, tn);
                    else if constexpr (std::is_same_v<ObjT, RenderTorus<T>>)
                        return intersectTorus(ray, o.shape, curMin, tMax, th, tn);
                    else if constexpr (std::is_same_v<ObjT, RenderQuadric<T>>)
                        return intersectQuadric(ray, o.shape, curMin, tMax, th, tn);
                    else if constexpr (std::is_same_v<ObjT, RenderSuperellipsoid<T>>)
                        return intersectSuperellipsoid(ray, o.shape, curMin, tMax, th, tn);
                    else if constexpr (std::is_same_v<ObjT, RenderDisc<T>>)
                        return intersectDisc(ray, o.shape, curMin, tMax, th, tn);
                    else if constexpr (std::is_same_v<ObjT, RenderHeightField<T>>)
                        return o.shape ? intersectHeightField(ray, *o.shape, curMin, tMax, th, tn) : false;
                    else if constexpr (std::is_same_v<ObjT, RenderBlob<T>>)
                        return o.shape ? intersectBlob(ray, *o.shape, curMin, tMax, th, tn) : false;
                    else if constexpr (std::is_same_v<ObjT, RenderLathe<T>>)
                        return o.shape ? intersectLathe(ray, *o.shape, curMin, tMax, th, tn) : false;
                    else if constexpr (std::is_same_v<ObjT, RenderPrism<T>>)
                        return o.shape ? intersectPrism(ray, *o.shape, curMin, tMax, th, tn) : false;
                    else if constexpr (std::is_same_v<ObjT, RenderMesh<T>>)
                        return o.shape ? intersectMesh(ray, *o.shape, curMin, tMax, th, tn) : false;
                    else if constexpr (std::is_same_v<ObjT, std::shared_ptr<CsgNode<T>>>)
                        return o ? intersectCsg(ray, *o, curMin, tMax, th, tn) : false;
                    else
                        return false;
                }, node.children[ci]);
                if (didHit && th < bestT) {
                    bestT = th; bestN = tn; bestChild = ci;
                }
            }
            if (bestChild < 0) return false;

            auto hitPoint = ray.at(bestT);
            bool valid = true;

            if (node.op == CsgNode<T>::Op::Intersection) {
                // 交差点が全子要素の内側にあるか
                for (int ci = 0; ci < static_cast<int>(node.children.size()); ++ci) {
                    if (ci == bestChild) continue;
                    if (!isInsideObject(hitPoint, node.children[ci])) {
                        valid = false; break;
                    }
                }
            } else { // Difference
                // 交差点が children[0] の内側かつ他の外側
                if (bestChild == 0) {
                    for (std::size_t ci = 1; ci < node.children.size(); ++ci) {
                        if (isInsideObject(hitPoint, node.children[ci])) {
                            valid = false; break;
                        }
                    }
                } else {
                    if (!isInsideObject(hitPoint, node.children[0])) valid = false;
                    for (std::size_t ci = 1; ci < node.children.size(); ++ci) {
                        if (static_cast<int>(ci) == bestChild) continue;
                        if (isInsideObject(hitPoint, node.children[ci])) {
                            valid = false; break;
                        }
                    }
                    if (valid) bestN = bestN * T(-1);  // 法線反転 (差の内面)
                }
            }

            if (valid) {
                tHit = bestT;
                normal = bestN;
                return true;
            }
            curMin = bestT + T(1e-4);  // 次の交差を探す
        }
        return false;
    }

} // namespace detail

// ========================================================================
// レイトレーサー本体
// ========================================================================

template<typename T>
class RayTracer {
public:
    RayTracer(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(60)) {
        camPos_ = position;
        camLookAt_ = lookAt;
        fov_ = fovDeg;
    }

    // --- シーン構築 ---
    void addLight(const Light<T>& light) {
        lights_.push_back(light);
    }

    /// 拡張光源を追加 (スポット, エリア, 平行光源)
    void addSceneLight(const SceneLight<T>& light) {
        sceneLights_.push_back(light);
    }

    /// スポットライト追加ヘルパー
    void addSpotLight(const Point3D<T>& pos, const Point3D<T>& dir,
                      const RGB<T>& color = {T(1),T(1),T(1)},
                      T falloff = T(45), T hotspot = T(30), T tightness = T(1)) {
        SceneLight<T> sl;
        sl.type = SceneLight<T>::Type::Spot;
        sl.position = pos; sl.direction = dir; sl.color = color;
        sl.falloffAngle = falloff; sl.hotspotAngle = hotspot; sl.tightness = tightness;
        sceneLights_.push_back(sl);
    }

    /// エリアライト追加ヘルパー
    void addAreaLight(const Point3D<T>& pos, const Point3D<T>& ax1,
                      const Point3D<T>& ax2, const RGB<T>& color = {T(1),T(1),T(1)},
                      int s1 = 4, int s2 = 4) {
        SceneLight<T> sl;
        sl.type = SceneLight<T>::Type::Area;
        sl.position = pos; sl.axis1 = ax1; sl.axis2 = ax2;
        sl.color = color; sl.samples1 = s1; sl.samples2 = s2;
        sceneLights_.push_back(sl);
    }

    /// 平行光源追加ヘルパー
    void addDirectionalLight(const Point3D<T>& dir,
                             const RGB<T>& color = {T(1),T(1),T(1)}) {
        SceneLight<T> sl;
        sl.type = SceneLight<T>::Type::Directional;
        sl.direction = dir; sl.color = color;
        sceneLights_.push_back(sl);
    }

    void addSphere(const Sphere<T>& s, const Material<T>& mat = {}) {
        objects_.push_back(RenderSphere<T>{s, mat});
    }

    void addPlane(const Plane<T>& p, const Material<T>& mat = {}) {
        objects_.push_back(RenderPlane<T>{p, mat});
    }

    void addTriangle(const Triangle3D<T>& tri, const Material<T>& mat = {}) {
        objects_.push_back(RenderTriangle<T>{tri, mat});
    }

    void addBox(const AABB<T>& box, const Material<T>& mat = {}) {
        objects_.push_back(RenderBox<T>{box, mat});
    }

    // --- POV-1: トーラス ---
    void addTorus(const Torus<T>& t, const Material<T>& mat = {}) {
        objects_.push_back(RenderTorus<T>{t, mat});
    }

    // --- POV-2: 二次曲面 ---
    void addQuadric(const Quadric<T>& q, const Material<T>& mat = {}) {
        objects_.push_back(RenderQuadric<T>{q, mat});
    }

    // --- POV-3: 超二次曲面 ---
    void addSuperellipsoid(const Superellipsoid<T>& se, const Material<T>& mat = {}) {
        objects_.push_back(RenderSuperellipsoid<T>{se, mat});
    }

    // --- POV-9: ディスク ---
    void addDisc(const Disc<T>& d, const Material<T>& mat = {}) {
        objects_.push_back(RenderDisc<T>{d, mat});
    }

    // --- POV-5: ハイトフィールド ---
    void addHeightField(const std::vector<std::vector<T>>& heights,
                        const AABB<T>& bounds, const Material<T>& mat = {}) {
        auto data = std::make_shared<HeightFieldData<T>>();
        data->heights = heights;
        data->bounds = bounds;
        objects_.push_back(RenderHeightField<T>{data, mat});
    }

    // --- POV-6: ブロブ ---
    void addBlob(const std::vector<BlobComponent<T>>& components,
                 T threshold = T(1), const Material<T>& mat = {}) {
        auto data = std::make_shared<BlobData<T>>();
        data->components = components;
        data->threshold = threshold;
        objects_.push_back(RenderBlob<T>{data, mat});
    }

    // --- POV-7: ラス (回転体) ---
    void addLathe(const std::vector<Point2D<T>>& points,
                  const Material<T>& mat = {}) {
        auto data = std::make_shared<LatheData<T>>();
        data->points = points;
        objects_.push_back(RenderLathe<T>{data, mat});
    }

    // --- POV-8: プリズム ---
    void addPrism(const std::vector<Point2D<T>>& points,
                  T y0, T y1, const Material<T>& mat = {}) {
        auto data = std::make_shared<PrismData<T>>();
        data->points = points;
        data->y0 = y0;
        data->y1 = y1;
        objects_.push_back(RenderPrism<T>{data, mat});
    }

    // --- POV-10: メッシュ ---
    void addMesh(const std::vector<Point3D<T>>& vertices,
                 const std::vector<std::array<int, 3>>& faces,
                 const std::vector<Point3D<T>>& normals = {},
                 const Material<T>& mat = {}) {
        auto data = std::make_shared<MeshData<T>>();
        data->vertices = vertices;
        data->faces = faces;
        data->normals = normals;
        data->computeBounds();
        objects_.push_back(RenderMesh<T>{data, mat});
    }

    // --- POV-4: CSG ---
    void addCsg(typename CsgNode<T>::Op op,
                const std::vector<SceneObject<T>>& children,
                const Material<T>& mat = {}) {
        auto node = std::make_shared<CsgNode<T>>();
        node->op = op;
        node->children = children;
        node->material = mat;
        objects_.push_back(node);
    }

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

    void setMaxDepth(int d) { maxDepth_ = d; }
    void setSamples(int n) { samples_ = n; }  // n×n 初期サンプリング
    void setAAThreshold(T t) { aaThreshold_ = t; }  // 適応的 AA 閾値
    void setMaxAADepth(int d) { maxAADepth_ = d; }   // 適応的再帰深度

    // --- レンダリング ---
    void render() {
        // カメラ座標系構築
        forward_ = (camLookAt_ - camPos_).normalized();
        right_ = cross(forward_, Point3D<T>{T(0), T(1), T(0)}).normalized();
        up_ = cross(right_, forward_);

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

        #pragma omp parallel for schedule(dynamic, 4)
        for (int py = 0; py < static_cast<int>(height_); ++py) {
            for (std::size_t px = 0; px < width_; ++px) {
                uint32_t rng = pixelSeed(px, static_cast<std::size_t>(py));
                auto color = adaptiveSample(
                    T(px), T(py), T(1), 0, 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;
            }
        }
    }

    const Image<T>& getImage() const { return image_; }
    Image<T>& getImage() { return image_; }

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

    Point3D<T> camPos_{T(0), T(5), T(-10)};
    Point3D<T> camLookAt_{T(0), T(0), T(0)};
    T fov_ = T(60);

    RGB<T> bgTop_{T(0.5), T(0.7), T(1.0)};
    RGB<T> bgBottom_{T(1.0), T(1.0), T(1.0)};

    int maxDepth_ = 5;
    int samples_ = 2;       // 初期グリッド n×n
    T aaThreshold_ = T(0.03); // 色差閾値 (小さいほど高品質)
    int maxAADepth_ = 3;     // 適応的再帰の最大深度
    std::vector<Light<T>> lights_;
    std::vector<SceneLight<T>> sceneLights_;
    std::vector<SceneObject<T>> objects_;

    // レンダリング時の一時変数 (render() で設定)
    Point3D<T> forward_, right_, up_;
    T halfW_, halfH_;

    // --- ピクセル座標ベースのハッシュ RNG (スレッドセーフ) ---
    static T pixelRand(uint32_t& state) {
        // xorshift32
        state ^= state << 13;
        state ^= state >> 17;
        state ^= state << 5;
        return T(state & 0xFFFFu) / T(0x10000u);  // [0, 1)
    }

    static uint32_t pixelSeed(std::size_t px, std::size_t py) {
        // ピクセル座標からシード生成 (hash)
        uint32_t h = uint32_t(px * 73856093u ^ py * 19349663u);
        h ^= h >> 16; h *= 0x45d9f3bu; h ^= h >> 16;
        return h | 1u;  // ゼロ回避
    }

    // --- ピクセル座標 → レイ → 色 ---
    RGB<T> sampleAt(T x, T y) 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();
        Ray<T> ray{camPos_, dir};
        return trace(ray, 0);
    }

    // --- RGB 差 (最大チャネル差) ---
    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)});
    }

    // --- 適応的スーパーサンプリング ---
    // (x, y) を左上とする size×size ピクセル領域を再帰的にサンプリング
    RGB<T> adaptiveSample(T x, T y, T size, int depth, uint32_t& rng) const {
        // 四隅 + ジッタサンプル
        T jit = size * T(0.1);  // ジッタ量
        RGB<T> c00 = sampleAt(x + jit * pixelRand(rng),
                              y + jit * pixelRand(rng));
        RGB<T> c10 = sampleAt(x + size + jit * (pixelRand(rng) - T(0.5)),
                              y + jit * pixelRand(rng));
        RGB<T> c01 = sampleAt(x + jit * pixelRand(rng),
                              y + size + jit * (pixelRand(rng) - T(0.5)));
        RGB<T> c11 = sampleAt(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 = sampleAt(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_) {
            // 4分割して再帰
            RGB<T> q00 = adaptiveSample(x,        y,        half, depth + 1, rng);
            RGB<T> q10 = adaptiveSample(x + half, y,        half, depth + 1, rng);
            RGB<T> q01 = adaptiveSample(x,        y + half, half, depth + 1, rng);
            RGB<T> q11 = adaptiveSample(x + half, y + half, half, depth + 1, rng);
            return (q00 + q10 + q01 + q11) * T(0.25);
        }

        // 5サンプルの平均
        return (c00 + c10 + c01 + c11 + cMid) * T(0.2);
    }

    // --- 背景グラデーション ---
    RGB<T> background(const Ray<T>& ray) const {
        T t = T(0.5) * (ray.direction.y + T(1));
        return (T(1) - t) * bgBottom_ + t * bgTop_;
    }

    // --- 最近交差を求める ---
    bool closestHit(const Ray<T>& ray, T tMin, T tMax,
                    HitRecord<T>& rec) const {
        bool hit = false;
        T closest = tMax;

        for (const auto& obj : objects_) {
            T tHit;
            Point3D<T> n;
            const Material<T>* mat = nullptr;
            bool didHit = false;

            std::visit([&](const auto& o) {
                using ObjT = std::decay_t<decltype(o)>;
                if constexpr (std::is_same_v<ObjT, std::shared_ptr<CsgNode<T>>>) {
                    if (o) {
                        mat = &o->material;
                        didHit = detail::intersectCsg(ray, *o, tMin, closest, tHit, n);
                    }
                } else {
                    mat = &o.material;
                    if constexpr (std::is_same_v<ObjT, RenderSphere<T>>)
                        didHit = detail::intersectSphere(ray, o.shape, tMin, closest, tHit, n);
                    else if constexpr (std::is_same_v<ObjT, RenderPlane<T>>)
                        didHit = detail::intersectPlane(ray, o.shape, tMin, closest, tHit, n);
                    else if constexpr (std::is_same_v<ObjT, RenderTriangle<T>>)
                        didHit = detail::intersectTriangle(ray, o.shape, tMin, closest, tHit, n);
                    else if constexpr (std::is_same_v<ObjT, RenderBox<T>>)
                        didHit = detail::intersectAABB(ray, o.shape, tMin, closest, tHit, n);
                    else if constexpr (std::is_same_v<ObjT, RenderTorus<T>>)
                        didHit = detail::intersectTorus(ray, o.shape, tMin, closest, tHit, n);
                    else if constexpr (std::is_same_v<ObjT, RenderQuadric<T>>)
                        didHit = detail::intersectQuadric(ray, o.shape, tMin, closest, tHit, n);
                    else if constexpr (std::is_same_v<ObjT, RenderSuperellipsoid<T>>)
                        didHit = detail::intersectSuperellipsoid(ray, o.shape, tMin, closest, tHit, n);
                    else if constexpr (std::is_same_v<ObjT, RenderDisc<T>>)
                        didHit = detail::intersectDisc(ray, o.shape, tMin, closest, tHit, n);
                    else if constexpr (std::is_same_v<ObjT, RenderHeightField<T>>)
                        didHit = o.shape ? detail::intersectHeightField(ray, *o.shape, tMin, closest, tHit, n) : false;
                    else if constexpr (std::is_same_v<ObjT, RenderBlob<T>>)
                        didHit = o.shape ? detail::intersectBlob(ray, *o.shape, tMin, closest, tHit, n) : false;
                    else if constexpr (std::is_same_v<ObjT, RenderLathe<T>>)
                        didHit = o.shape ? detail::intersectLathe(ray, *o.shape, tMin, closest, tHit, n) : false;
                    else if constexpr (std::is_same_v<ObjT, RenderPrism<T>>)
                        didHit = o.shape ? detail::intersectPrism(ray, *o.shape, tMin, closest, tHit, n) : false;
                    else if constexpr (std::is_same_v<ObjT, RenderMesh<T>>)
                        didHit = o.shape ? detail::intersectMesh(ray, *o.shape, tMin, closest, tHit, n) : false;
                }
            }, obj);

            if (didHit) {
                hit = true;
                closest = tHit;
                rec.t = tHit;
                rec.point = ray.at(tHit);
                rec.material = mat;
                rec.frontFace = dot(ray.direction, n) < T(0);
                rec.normal = rec.frontFace ? n : (T(-1) * n);
            }
        }
        return hit;
    }

    // --- シャドウ判定 ---
    bool inShadow(const Point3D<T>& point, const Point3D<T>& lightPos) const {
        auto toLight = lightPos - point;
        T dist = toLight.norm();
        auto dir = toLight * (T(1) / dist);
        Ray<T> shadowRay{point + T(1e-4) * dir, dir};
        HitRecord<T> tmp;
        return closestHit(shadowRay, T(1e-4), dist - T(1e-4), tmp);
    }

    // --- 反射方向 ---
    static Point3D<T> reflect(const Point3D<T>& v, const Point3D<T>& n) {
        return v - T(2) * dot(v, n) * n;
    }

    // --- 屈折方向 (Snell) ---
    static bool refract(const Point3D<T>& v, const Point3D<T>& n,
                        T etaRatio, Point3D<T>& refracted) {
        T cosI = -dot(v, n);
        T sin2T = etaRatio * etaRatio * (T(1) - cosI * cosI);
        if (sin2T > T(1)) return false;  // 全反射
        T cosT = std::sqrt(T(1) - sin2T);
        refracted = etaRatio * v + (etaRatio * cosI - cosT) * n;
        return true;
    }

    // --- Schlick の Fresnel 近似 ---
    static T schlick(T cosine, T ior) {
        T r0 = (T(1) - ior) / (T(1) + ior);
        r0 = r0 * r0;
        return r0 + (T(1) - r0) * std::pow(T(1) - cosine, T(5));
    }

    // --- Phong シェーディング (テクスチャ・バンプ対応) ---
    RGB<T> shade(const Ray<T>& ray, const HitRecord<T>& rec) const {
        const auto& mat = *rec.material;
        RGB<T> surfColor = mat.colorAt(rec.point);
        Point3D<T> N = mat.normalAt(rec.point, rec.normal);
        RGB<T> result = surfColor * mat.ambient;

        // 虹色干渉 (iridescence)
        if (mat.iridescence > T(0)) {
            T cosI = std::abs(dot(ray.direction, N));
            T hue = std::fmod(cosI * T(6), T(1));
            // 簡易 HSV→RGB (S=1, V=iridescence)
            T h6 = hue * T(6);
            int hi = static_cast<int>(h6) % 6;
            T f = h6 - std::floor(h6);
            T v = mat.iridescence;
            RGB<T> iridColor;
            switch (hi) {
                case 0: iridColor = {v, v*f, T(0)}; break;
                case 1: iridColor = {v*(T(1)-f), v, T(0)}; break;
                case 2: iridColor = {T(0), v, v*f}; break;
                case 3: iridColor = {T(0), v*(T(1)-f), v}; break;
                case 4: iridColor = {v*f, T(0), v}; break;
                default: iridColor = {v, T(0), v*(T(1)-f)}; break;
            }
            surfColor = surfColor + iridColor;
        }

        auto V = (T(-1) * ray.direction);

        // 1光源分のシェーディング計算
        auto shadeOneLight = [&](const Point3D<T>& L, const RGB<T>& lc, T atten) {
            T diff = std::max(T(0), dot(N, L));
            result += mulRgb(surfColor, lc) * (mat.diffuse * diff * atten);
            auto H = (L + V).normalized();
            T spec = std::pow(std::max(T(0), dot(N, H)), mat.shininess);
            if (mat.metallic)
                result += mulRgb(surfColor, lc) * (mat.specular * spec * atten);
            else
                result += lc * (mat.specular * spec * atten);
        };

        // 従来の点光源
        for (const auto& light : lights_) {
            if (inShadow(rec.point, light.position)) continue;
            auto L = (light.position - rec.point).normalized();
            shadeOneLight(L, light.color, T(1));
        }

        // 拡張光源
        for (const auto& sl : sceneLights_) {
            if (sl.type == SceneLight<T>::Type::Point) {
                if (inShadow(rec.point, sl.position)) continue;
                auto L = (sl.position - rec.point).normalized();
                shadeOneLight(L, sl.color, sl.intensity);
            }
            else if (sl.type == SceneLight<T>::Type::Spot) {
                if (inShadow(rec.point, sl.position)) continue;
                T atten = sl.spotAttenuation(rec.point) * sl.intensity;
                if (atten < T(1e-4)) continue;
                auto L = (sl.position - rec.point).normalized();
                shadeOneLight(L, sl.color, atten);
            }
            else if (sl.type == SceneLight<T>::Type::Directional) {
                auto L = (T(-1) * sl.direction).normalized();
                // 平行光源: 無限遠からのシャドウチェック
                Ray<T> shadowRay{rec.point + T(1e-4) * L, L};
                HitRecord<T> tmp;
                if (closestHit(shadowRay, T(1e-4), T(1e10), tmp)) continue;
                shadeOneLight(L, sl.color, sl.intensity);
            }
            else if (sl.type == SceneLight<T>::Type::Area) {
                // エリアライト: 複数サンプルの平均 (ソフトシャドウ)
                uint32_t rng = static_cast<uint32_t>(
                    rec.point.x * T(73856093) + rec.point.z * T(19349663));
                rng |= 1u;
                auto samples = sl.areaSamples(rng);
                T contrib = T(0);
                int lit = 0;
                for (const auto& sp : samples) {
                    if (!inShadow(rec.point, sp)) {
                        auto L = (sp - rec.point).normalized();
                        shadeOneLight(L, sl.color,
                                      sl.intensity / T(samples.size()));
                        ++lit;
                    }
                }
            }
        }

        return result;
    }

    // --- 再帰レイトレーシング ---
    RGB<T> trace(const Ray<T>& ray, int depth) const {
        if (depth >= maxDepth_) return RGB<T>{T(0), T(0), T(0)};

        HitRecord<T> rec;
        if (!closestHit(ray, T(1e-4), T(1e10), rec))
            return background(ray);

        const auto& mat = *rec.material;

        // Phong シェーディング
        RGB<T> color = shade(ray, rec);

        // 反射
        if (mat.reflectivity > T(0) && depth < maxDepth_) {
            auto reflDir = reflect(ray.direction, rec.normal).normalized();
            Ray<T> reflRay{rec.point + T(1e-4) * reflDir, reflDir};
            auto reflColor = trace(reflRay, depth + 1);
            color = color * (T(1) - mat.reflectivity) + reflColor * mat.reflectivity;
        }

        // 屈折
        if (mat.transparency > T(0) && depth < maxDepth_) {
            T etaRatio = rec.frontFace ? (T(1) / mat.ior) : mat.ior;
            Point3D<T> refracted;
            T cosI = std::min(T(1), std::abs(dot(ray.direction, rec.normal)));
            T fresnel = schlick(cosI, mat.ior);

            if (refract(ray.direction, rec.normal, etaRatio, refracted)) {
                auto refrDir = refracted.normalized();
                Ray<T> refrRay{rec.point + T(1e-4) * refrDir, refrDir};
                auto refrColor = trace(refrRay, depth + 1);
                // Fresnel による反射/屈折の混合
                color = color * fresnel + refrColor * (mat.transparency * (T(1) - fresnel));
            }
            // 全反射の場合は反射のみ (上の reflectivity で処理済み)
        }

        return color;
    }

    // RGB の要素ごと乗算ヘルパー
    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};
    }
};

} // namespace calx::image

#endif // CALX_IMAGE_RENDER_RAYTRACER_HPP
