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

/**
 * @file mandelbulb.hpp
 * @brief Mandelbulb フラクタル レイマーチング レンダラー (DRW-70)
 *
 * 3D Mandelbrot 集合の一般化である Mandelbulb を
 * distance estimation + レイマーチングで描画する。
 *
 * 使用例:
 *   using namespace calx::image;
 *   MandelbulbRenderer<double> mb(800, 600);
 *   mb.camera({0, 0, -3}, {0, 0, 0});
 *   mb.render();
 *   saveBmp("mandelbulb.bmp", mb.image());
 */

#ifndef CALX_IMAGE_RENDER_MANDELBULB_HPP
#define CALX_IMAGE_RENDER_MANDELBULB_HPP

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

namespace calx::image {

/// Mandelbulb レンダラー
template<typename T = double>
class MandelbulbRenderer {
public:
    /// コンストラクタ
    MandelbulbRenderer(int width, int height)
        : width_(width), height_(height), img_(height, width, 3) {}

    /// カメラ設定
    void camera(const Point3D<T>& eye, const Point3D<T>& target,
                const Point3D<T>& up = {T(0), T(1), T(0)}, T fovDeg = T(60)) {
        eye_    = eye;
        target_ = target;
        up_     = up;
        fovDeg_ = fovDeg;
    }

    /// Mandelbulb パラメータ
    void setParams(T power = T(8), int maxIter = 12, T bailout = T(2)) {
        power_   = power;
        maxIter_ = maxIter;
        bailout_ = bailout;
    }

    /// レイマーチング パラメータ
    void setRayParams(int maxSteps = 256, T minDist = T(1e-4), T maxDist = T(10)) {
        maxSteps_ = maxSteps;
        minDist_  = minDist;
        maxDist_  = maxDist;
    }

    /// 光源方向
    void setLightDir(const Point3D<T>& dir) { lightDir_ = dir.normalized(); }

    /// 背景色
    void setBackground(const RGB<T>& bg) { bg_ = bg; }

    /// スーパーサンプリング倍率 (1=なし, 2=2x2, 3=3x3)
    void setSuperSample(int n) { superSample_ = std::max(1, n); }

    /// レンダリング実行
    void render() {
        auto forward = (target_ - eye_).normalized();
        auto right   = cross(forward, up_).normalized();
        auto camUp   = cross(right, forward);
        T halfH      = std::tan(fovDeg_ * T(3.14159265358979323846) / T(360));
        T aspect     = T(width_) / T(height_);
        int ss       = superSample_;
        T ssInv      = T(1) / T(ss * ss);

        #pragma omp parallel for schedule(dynamic, 4)
        for (int y = 0; y < height_; ++y) {
            for (int x = 0; x < width_; ++x) {
                RGB<T> accum{T(0), T(0), T(0)};
                for (int sy = 0; sy < ss; ++sy) {
                    for (int sx = 0; sx < ss; ++sx) {
                        T px = T(x) + (T(sx) + T(0.5)) / T(ss);
                        T py = T(y) + (T(sy) + T(0.5)) / T(ss);
                        T u = (T(2) * px / T(width_)  - T(1)) * halfH * aspect;
                        T v = (T(1) - T(2) * py / T(height_)) * halfH;
                        auto dir = (forward + right * u + camUp * v).normalized();
                        accum += trace(eye_, dir);
                    }
                }
                auto color = accum * ssInv;
                img_.at(y, x, 0) = color.r;
                img_.at(y, x, 1) = color.g;
                img_.at(y, x, 2) = color.b;
            }
        }
    }

    /// 結果画像を取得
    const Image<T>& image() const { return img_; }

private:
    int width_, height_;
    Image<T> img_;

    // カメラ
    Point3D<T> eye_{T(0), T(0), T(-3)};
    Point3D<T> target_{T(0), T(0), T(0)};
    Point3D<T> up_{T(0), T(1), T(0)};
    T fovDeg_ = T(60);

    // Mandelbulb
    T power_   = T(8);
    int maxIter_ = 12;
    T bailout_ = T(2);

    // レイマーチング
    int maxSteps_ = 256;
    T minDist_  = T(1e-4);
    T maxDist_  = T(10);

    // スーパーサンプリング
    int superSample_ = 1;

    // ライティング
    Point3D<T> lightDir_ = Point3D<T>(T(-0.5), T(1), T(-0.3)).normalized();
    RGB<T> bg_{T(0.05), T(0.05), T(0.1)};

    /// Mandelbulb の距離推定関数 (DE)
    /// @return {距離, 反復回数の比率 [0,1]}
    struct DEResult { T dist; T iterRatio; };

    DEResult distanceEstimate(const Point3D<T>& pos) const {
        Point3D<T> z = pos;
        T dr = T(1);    // ∂|z|/∂c の追跡
        T r  = T(0);
        int i = 0;

        for (; i < maxIter_; ++i) {
            r = z.norm();
            if (r > bailout_) break;

            // 球面座標変換
            T theta = std::acos(std::clamp(z.z / r, T(-1), T(1)));
            T phi   = std::atan2(z.y, z.x);
            T rn    = std::pow(r, power_);

            // ∂r^n/∂c の更新
            dr = std::pow(r, power_ - T(1)) * power_ * dr + T(1);

            // z = z^power + c
            T sinT = std::sin(power_ * theta);
            T cosT = std::cos(power_ * theta);
            T sinP = std::sin(power_ * phi);
            T cosP = std::cos(power_ * phi);

            z.x = pos.x + rn * sinT * cosP;
            z.y = pos.y + rn * sinT * sinP;
            z.z = pos.z + rn * cosT;
        }

        T dist = T(0.5) * std::log(r) * r / dr;

        // 連続的着色 (smooth iteration count)
        // i + 1 - log2(log(|z|) / log(bailout)) で反復回数を補間
        T smooth = T(i);
        if (i < maxIter_ && r > T(1)) {
            smooth = T(i) + T(1) - std::log(std::log(r) / std::log(bailout_))
                                  / std::log(power_);
        }
        T iterRatio = std::clamp(smooth / T(maxIter_), T(0), T(1));

        return {dist, iterRatio};
    }

    /// 法線推定 (中心差分)
    Point3D<T> estimateNormal(const Point3D<T>& p) const {
        T eps = T(1e-5);
        T dx = distanceEstimate({p.x + eps, p.y, p.z}).dist
             - distanceEstimate({p.x - eps, p.y, p.z}).dist;
        T dy = distanceEstimate({p.x, p.y + eps, p.z}).dist
             - distanceEstimate({p.x, p.y - eps, p.z}).dist;
        T dz = distanceEstimate({p.x, p.y, p.z + eps}).dist
             - distanceEstimate({p.x, p.y, p.z - eps}).dist;
        return Point3D<T>(dx, dy, dz).normalized();
    }

    /// 法線ベースのマテリアル色 (滑らか、ノイズなし)
    RGB<T> materialColor(const Point3D<T>& normal) const {
        // 法線方向で色を決定 (連続的で滑らか)
        // 上向き: 暖色 (琥珀), 横向き: 青銅, 下向き: 深い藍
        T ny = normal.y * T(0.5) + T(0.5);  // [0, 1]
        T nx = normal.x * T(0.5) + T(0.5);

        RGB<T> warm{T(0.85), T(0.55), T(0.25)};   // 琥珀
        RGB<T> mid{T(0.45), T(0.35), T(0.55)};     // 青銅紫
        RGB<T> cool{T(0.15), T(0.20), T(0.45)};    // 深い藍

        // 上下で暖色⇔寒色
        auto base = warm * ny + cool * (T(1) - ny);
        // 左右で紫味を加える
        base = base * (T(1) - nx * T(0.2)) + mid * (nx * T(0.2));
        return base;
    }

    /// 距離に応じたフォグ
    RGB<T> applyFog(const RGB<T>& color, T totalDist) const {
        T fogAmount = T(1) - std::exp(-totalDist * T(0.05));
        return color * (T(1) - fogAmount) + bg_ * fogAmount;
    }

    /// レイマーチング + シェーディング
    RGB<T> trace(const Point3D<T>& origin, const Point3D<T>& dir) const {
        T totalDist = T(0);
        T minDE = T(1e10);

        for (int step = 0; step < maxSteps_; ++step) {
            auto p = origin + dir * totalDist;
            auto [d, iterRatio] = distanceEstimate(p);
            minDE = std::min(minDE, d);

            if (d < minDist_) {
                auto normal = estimateNormal(p);

                // --- 4 光源 ---
                T ndotl1 = std::max(T(0), dot(normal, lightDir_));
                auto light2Dir = Point3D<T>(T(0.3), T(-0.5), T(0.8)).normalized();
                T ndotl2 = std::max(T(0), dot(normal, light2Dir));
                auto light3Dir = Point3D<T>(T(-0.8), T(0.2), T(0.5)).normalized();
                T ndotl3 = std::max(T(0), dot(normal, light3Dir));
                // カメラ光源 (ヘッドライト): レイの逆方向
                T ndotCam = std::max(T(0), dot(normal, T(-1) * dir));

                // AO (適応的ステップサイズ)
                T ao = T(1);
                T aoScale = std::max(minDist_ * T(10), T(0.005));
                for (int k = 1; k <= 5; ++k) {
                    T dist = aoScale * T(k);
                    auto aoP = p + normal * dist;
                    T aoDist = distanceEstimate(aoP).dist;
                    ao -= std::max(T(0), dist - aoDist) / (dist * T(5));
                }
                ao = std::clamp(ao, T(0.15), T(1));

                // マテリアル色 (法線ベース = 滑らか)
                auto baseColor = materialColor(normal);

                // アンビエント (半球ライティング、強め)
                T ambient = T(0.30) + T(0.10) * normal.y;

                // ディフューズ 4 光源
                RGB<T> keyLight  = baseColor * (ndotl1 * T(1.0));
                RGB<T> fillLight = baseColor * RGB<T>(T(0.3), T(0.4), T(0.7)) * (ndotl2 * T(0.4));
                RGB<T> backLight = baseColor * RGB<T>(T(0.6), T(0.5), T(0.3)) * (ndotl3 * T(0.3));
                RGB<T> camLight  = baseColor * RGB<T>(T(0.9), T(0.9), T(1.0)) * (ndotCam * T(0.5));

                // スペキュラ
                auto halfVec = (lightDir_ + (T(-1) * dir)).normalized();
                T ndoth = std::max(T(0), dot(normal, halfVec));
                T spec = T(0.4) * std::pow(ndoth, T(40));

                // カメラ光源スペキュラ (常に正面から)
                T specCam = T(0.2) * std::pow(ndotCam, T(24));

                // フレネル リムライト
                T fresnel = std::pow(T(1) - ndotCam, T(4));
                T rim = fresnel * T(0.4);

                RGB<T> result = baseColor * (ambient * ao)
                              + (keyLight + fillLight + backLight + camLight) * ao
                              + RGB<T>(T(1), T(0.95), T(0.9)) * (spec + specCam)
                              + RGB<T>(T(0.4), T(0.5), T(0.8)) * rim;

                return applyFog(result, totalDist).clamped();
            }

            totalDist += d;
            if (totalDist > maxDist_) break;
        }

        T glow = std::clamp(T(1) / (T(1) + minDE * T(500)), T(0), T(1));
        glow = glow * glow;
        RGB<T> glowColor{T(0.1) * glow, T(0.12) * glow, T(0.25) * glow};
        return (bg_ + glowColor).clamped();
    }
};

} // namespace calx::image

#endif // CALX_IMAGE_RENDER_MANDELBULB_HPP
