#pragma once
/**
 * @file sdf_optimize.h
 * @brief 基于梯度的 SDF 形状优化
 *
 * 提供 SDF 形状的拟合、碰撞解析、可及性分析和对称性检测等功能。
 *
 * ## 主要功能
 *
 * | 功能                   | 说明                                             |
 * |------------------------|-------------------------------------------------|
 * | Shape Fitting          | 将参数化形状拟合到点云或目标 SDF                  |
 * | Collision Avoidance    | 平移形状以解决穿透                                |
 * | Accessibility          | 采样半球以评估表面上点的可及性                     |
 * | Symmetry Detection     | 检测反射对称性并找到最佳对称平面                   |
 * | Volume & Mass          | 通过蒙地卡罗法估计体积和质心                       |
 * | Parameter Optimization | 收集所有可优化的参数并计算数值梯度                  |
 *
 * ## 使用模式
 *
 * ```cpp
 * auto shape = SdfNode::sphere(0.5);
 * auto result = fit_to_point_cloud(shape, points, 0.01, 200);
 * if (result.converged) {
 *     // result.optimized_shape 包含优化后的参数
 * }
 * ```
 *
 * @ingroup sdf
 */
#include "vde/sdf/sdf_tree.h"
#include "vde/core/point.h"
#include <vector>
#include <functional>

namespace vde::sdf {

using core::Point3D;
using core::Vector3D;

// ────────────────────────────────────────────────
//  Shape Optimization
// ────────────────────────────────────────────────

/**
 * @brief 优化结果
 *
 * 包含优化后的形状、损失值、收敛标志和迭代历史。
 */
struct OptimizeResult {
    SdfNodePtr optimized_shape;      ///< 优化后的 SDF 树（参数已更新）
    double final_loss;               ///< 最终损失值
    int iterations;                  ///< 实际迭代次数
    bool converged;                  ///< 是否收敛
    std::vector<double> loss_history; ///< 每次迭代的损失记录
};

/**
 * @brief 损失函数类型
 *
 * 接受一个 SDF 求值函数（封装了参数化的形状），返回标量损失。
 *
 * @see fit_to_point_cloud 使用该类型的内置实现
 */
using LossFn = std::function<double(const std::function<double(const Point3D&)>&)>;

// ── Shape Fitting ───────────────────────────────

/**
 * @brief 将参数化形状拟合到目标点云
 *
 * 最小化目标表面点到拟合形状表面的有符号距离。
 *
 * **算法概要:**
 * 1. 遍历 target_points 计算每个点到初始形状的 SDF 值
 * 2. 用有限差分计算形状参数对损失的梯度
 * 3. 梯度下降更新参数
 * 4. 重复直至收敛或达到最大迭代次数
 *
 * @param initial 初始形状（sphere, box, cylinder 等），其参数将被优化
 * @param target_points 目标点云（表面采样点）
 * @param learning_rate 梯度下降步长（默认 0.01）
 * @param max_iterations 最大迭代次数（默认 100）
 * @return 优化结果（含优化后的形状和收敛信息）
 *
 * @note 损失函数为 L = mean(|sdf(p_i)|)，即点到表面的平均绝对距离。
 * @note initial 的共享子节点不会被复制——返回的 optimize_shape 共享未修改的子树。
 *
 * @code{.cpp}
 * std::vector<Point3D> points = { point cloud data };
 * auto sphere = SdfNode::sphere(1.0);
 * auto result = fit_to_point_cloud(sphere, points, 0.01, 200);
 * // result.optimized_shape 的 params.radius 已被优化
 * @endcode
 *
 * @see fit_surface_to_points 类似但仅最小化 |sdf| 而非有符号距离
 * @see fit_to_sdf 拟合到另一个 SDF
 */
[[nodiscard]] OptimizeResult fit_to_point_cloud(
    const SdfNodePtr& initial,
    const std::vector<Point3D>& target_points,
    double learning_rate = 0.01,
    int max_iterations = 100);

/**
 * @brief 将形状表面移动到目标采样点
 *
 * 与 fit_to_point_cloud 类似，但损失函数使用 L = mean(|sdf(p_i)|)。
 * 这强制采样点位于表面上（而非仅靠近表面）。
 *
 * @param initial 初始形状
 * @param surface_points 目标表面点
 * @param learning_rate 学习率
 * @param max_iterations 最大迭代次数
 * @return 优化结果
 *
 * @see fit_to_point_cloud 更通用的点云拟合（可包含内部/外部点）
 */
[[nodiscard]] OptimizeResult fit_surface_to_points(
    const SdfNodePtr& initial,
    const std::vector<Point3D>& surface_points,
    double learning_rate = 0.01,
    int max_iterations = 100);

/**
 * @brief 拟合一个 SDF 形状以匹配另一个目标 SDF
 *
 * 在规则网格上采样两个 SDF，最小化它们之间的均方差。
 *
 * **损失函数:** L = mean((sdf_source(p_i) - sdf_target(p_i))²)
 *
 * 适合将一个参数化形状"画像"成一个更复杂的隐式形状。
 *
 * @param source 被优化的参数化形状
 * @param target_sdf 目标 SDF 函数
 * @param bmin 采样网格最小角点
 * @param bmax 采样网格最大角点
 * @param grid_resolution 采样网格分辨率（默认 16）
 * @param learning_rate 学习率
 * @param max_iterations 最大迭代次数
 * @return 优化结果
 *
 * @warning 网格分辨率每增加一倍，计算量增加 8 倍。使用 8~32 可获得合理速度。
 *
 * @see fit_to_point_cloud 拟合到点云
 */
[[nodiscard]] OptimizeResult fit_to_sdf(
    const SdfNodePtr& source,
    const std::function<double(const Point3D&)>& target_sdf,
    const Point3D& bmin, const Point3D& bmax,
    int grid_resolution = 16,
    double learning_rate = 0.01,
    int max_iterations = 100);

// ── Collision Avoidance ─────────────────────────

/**
 * @brief 通过平移解决两个形状之间的穿透
 *
 * 在穿透区域采样点，沿梯度方向平移形状以分离它们。
 * 修改传入的形状（原地修改 translate_offset 参数）。
 *
 * @param shape_a 形状 A（原地修改）
 * @param shape_b 形状 B（原地修改）
 * @param step_size 每次迭代的平移步长（默认 0.1）
 * @param max_iterations 最大迭代次数（默认 50）
 * @return true 如果碰撞被成功解决；false 如果在迭代次数内未能解决
 *
 * @code{.cpp}
 * auto sphere_a = SdfNode::sphere(1.0);
 * auto sphere_b = SdfNode::translate(SdfNode::sphere(1.0), Point3D(1.5, 0, 0));
 * bool resolved = resolve_collision(sphere_a, sphere_b, 0.1, 100);
 * @endcode
 *
 * @see penetration_depth 估算穿透深度
 */
[[nodiscard]] bool resolve_collision(
    SdfNodePtr& shape_a, SdfNodePtr& shape_b,
    double step_size = 0.1,
    int max_iterations = 50);

/**
 * @brief 估算两个形状之间的最小穿透距离
 *
 * 在包围盒内随机采样，找到两个 SDF 均 ≤ 0 的区域（重叠区），
 * 计算重叠点到表面的最大距离。
 *
 * @param shape_a 形状 A
 * @param shape_b 形状 B
 * @param bmin 公共包围盒最小角点
 * @param bmax 公共包围盒最大角点
 * @param samples 蒙地卡罗采样数（默认 1000）
 * @return 最大穿透深度（负值表示不重叠时的最小间距）
 *
 * @see resolve_collision 解决穿透
 */
[[nodiscard]] double penetration_depth(
    const SdfNodePtr& shape_a,
    const SdfNodePtr& shape_b,
    const Point3D& bmin, const Point3D& bmax,
    int samples = 1000);

// ── Reachability / Accessibility ─────────────────

/**
 * @brief 计算形状表面上点的可及性分数
 *
 * 从表面点 p 发射半球光线，统计未被形状遮挡的比例。
 * 返回 0（完全不可及）到 1（完全可及）之间的分数。
 *
 * @param shape SDF 形状
 * @param p 表面点（应满足 |sdf(p)| < ε）
 * @param direction 接近方向
 * @param hemisphere_samples 半球采样光线数量（默认 64）
 * @return 可及性分数 [0, 1]
 *
 * @note 采样数越多越精确，但线性增加计算量。
 *
 * @see find_accessible_point 找到最大可及性的表面点
 */
[[nodiscard]] double accessibility(
    const SdfNodePtr& shape,
    const Point3D& p,           // surface point
    const Vector3D& direction,  // approach direction
    int hemisphere_samples = 64);

/**
 * @brief 找到给定接近方向下具有最大可及性的表面点
 *
 * 在包围盒内的规则网格上搜索，计算每个网格点处的可及性分数。
 *
 * @param shape SDF 形状
 * @param approach_dir 接近方向（如工具接近方向）
 * @param bmin 搜索包围盒最小角点
 * @param bmax 搜索包围盒最大角点
 * @param grid_res 搜索网格分辨率（默认 32）
 * @return 最大可及性的表面点
 *
 * @see accessibility 单点可及性计算
 */
[[nodiscard]] Point3D find_accessible_point(
    const SdfNodePtr& shape,
    const Vector3D& approach_dir,
    const Point3D& bmin, const Point3D& bmax,
    int grid_res = 32);

// ── Symmetry Detection ──────────────────────────

/**
 * @brief 检测形状在给定方向上的反射对称性
 *
 * 在形状的包围盒内随机采样点对 (p, p') 其中 p' 是关于给定平面的反射。
 * 比较两点的 SDF 值以评估对称程度。
 *
 * @param shape SDF 形状
 * @param bmin 包围盒最小角点
 * @param bmax 包围盒最大角点
 * @param plane_normal 反射平面的法向量
 * @param plane_offset 反射平面的偏移量（沿法向量方向，默认 0）
 * @param samples 采样点对数（默认 1000）
 * @return 对称性分数 [0, 1]（0 = 完全不对称，1 = 完美对称）
 *
 * @see find_symmetry_plane 自动寻找最佳对称平面
 */
[[nodiscard]] double symmetry_score(
    const SdfNodePtr& shape,
    const Point3D& bmin, const Point3D& bmax,
    const Vector3D& plane_normal,
    double plane_offset = 0.0,
    int samples = 1000);

/**
 * @brief 对称检测结果
 */
struct SymmetryResult {
    double score;      ///< 对称性分数 [0, 1]
    Vector3D normal;   ///< 最佳对称平面的法向量
    double offset;     ///< 最佳对称平面的偏移量
};

/**
 * @brief 自动找到形状的最佳反射对称平面
 *
 * 尝试多个候选平面方向（沿主轴和若干采样方向），
 * 返回对称性分数最高的平面。
 *
 * @param shape SDF 形状
 * @param bmin 包围盒最小角点
 * @param bmax 包围盒最大角点
 * @param samples 每个候选方向的采样数（默认 1000）
 * @return 最佳对称平面及其分数
 *
 * @see symmetry_score 给定方向的对称性检测
 */
[[nodiscard]] SymmetryResult find_symmetry_plane(
    const SdfNodePtr& shape,
    const Point3D& bmin, const Point3D& bmax,
    int samples = 1000);

// ── SDF Volume Computation ──────────────────────

/**
 * @brief 通过蒙地卡罗法估算隐式形状的体积
 *
 * 在包围盒内均匀随机采样，统计落在形状内部的点比例，
 * 乘以包围盒体积得到近似体积。
 *
 * **公式:** V ≈ V_box · N_inside / N_total，其中 N_inside = count(sdf(p_i) < 0)
 *
 * @param shape SDF 形状
 * @param bmin 采样包围盒最小角点
 * @param bmax 采样包围盒最大角点
 * @param samples 蒙地卡罗采样数（默认 10000）
 * @return 近似体积
 *
 * @note 精度与 √samples 成正比。10000 采样 ≈ 1% 相对误差。
 * @note 包围盒应尽可能紧致以减小方差。
 *
 * @see center_of_mass 蒙地卡罗质心估算
 */
[[nodiscard]] double estimate_volume(
    const SdfNodePtr& shape,
    const Point3D& bmin, const Point3D& bmax,
    int samples = 10000);

/**
 * @brief 通过蒙地卡罗法估算形状的质心
 *
 * 对形状内部的采样点取算术平均。
 *
 * **公式:** COM ≈ Σ p_i / N_inside，其中 sdf(p_i) < 0
 *
 * @param shape SDF 形状
 * @param bmin 采样包围盒最小角点
 * @param bmax 采样包围盒最大角点
 * @param samples 蒙地卡罗采样数（默认 10000）
 * @return 近似质心坐标
 *
 * @see estimate_volume 体积估算
 */
[[nodiscard]] Point3D center_of_mass(
    const SdfNodePtr& shape,
    const Point3D& bmin, const Point3D& bmax,
    int samples = 10000);

// ── Parameter Collection ────────────────────────

/**
 * @brief SDF 树中可变参数的引用
 *
 * 指向 SdfNode::params 中具体数值字段的指针，
 * 配合 collect_params() 收集所有可优化参数。
 */
struct ParamRef {
    double* value;     ///< 指向参数值的可写指针
    std::string name;  ///< 参数名称（调试用）
};

/**
 * @brief 从叶节点收集所有可变数值参数
 *
 * 遍历 SDF 树，收集所有图元节点中可优化的参数引用。
 * 返回的 ParamRef 列表可传递给 numerical_gradient() 进行优化。
 *
 * @param root SDF 树根节点（可写引用，因为参数将被原地修改）
 * @return 可变参数引用列表
 *
 * @code{.cpp}
 * auto root = SdfNode::op_union(
 *     SdfNode::sphere(1.0),
 *     SdfNode::box(Point3D(0.5, 0.5, 0.5))
 * );
 * auto params = collect_params(root);
 * // params[0].value 指向球体半径
 * // params[1..4] 指向盒子的 extents 分量
 * @endcode
 *
 * @see numerical_gradient 计算损失对参数的数值梯度
 */
[[nodiscard]] std::vector<ParamRef> collect_params(SdfNodePtr& root);

/**
 * @brief 计算标量损失对收集到的参数的数值梯度
 *
 * 对每个参数逐个做中心差分，返回梯度向量。
 *
 * @param params collect_params() 返回的参数引用列表
 * @param loss_fn 损失函数（无参数的可调用对象）
 * @param eps 差分步长（默认 1e-6）
 * @return 梯度向量 grad_loss，grad_loss[i] = ∂loss / ∂params[i].value
 *
 * @code{.cpp}
 * auto params = collect_params(root);
 * auto loss = [&]() {
 *     double s = 0;
 *     for (auto& pt : target_points) s += std::abs(evaluate(root, pt));
 *     return s / target_points.size();
 * };
 * auto grad = numerical_gradient(params, loss, 1e-6);
 * @endcode
 *
 * @see collect_params 收集参数
 * @see GradientDescent 梯度下降优化器
 */
[[nodiscard]] std::vector<double> numerical_gradient(
    const std::vector<ParamRef>& params,
    const std::function<double()>& loss_fn,
    double eps = 1e-6);

// ── Gradient Descent Utility ────────────────────

/**
 * @brief 固定学习率的简单梯度下降优化器
 *
 * 管理迭代计数和学习率，提供 step() 方法执行单步参数更新。
 *
 * @code{.cpp}
 * auto params = collect_params(root);
 * GradientDescent gd(0.01);
 * for (int i = 0; i < 100; ++i) {
 *     auto loss = gd.step(param_vals, [&](auto& vals, auto& grad) {
 *         // 填充 grad，返回 loss
 *     });
 *     if (loss < tol) break;
 * }
 * @endcode
 */
class GradientDescent {
public:
    /**
     * @brief 构造优化器
     * @param lr 学习率（步长因子）
     */
    explicit GradientDescent(double lr) : lr_(lr) {}

    /**
     * @brief 执行单步梯度下降
     *
     * 调用 grad_fn 获取当前参数的梯度和损失值，
     * 然后沿负梯度方向更新参数: param_i -= lr · grad_i。
     *
     * @param params 参数向量（原地修改）
     * @param grad_fn 梯度计算函数，签名为
     *        double(const std::vector<double>& vals, std::vector<double>& grad_out)
     *        接受当前参数值，填充 grad_out 并返回损失标量。
     * @return 当前损失值
     */
    double step(std::vector<double>& params,
                const std::function<double(const std::vector<double>&,
                                           std::vector<double>&)>& grad_fn);

    /** @brief 设置学习率 */
    void set_learning_rate(double lr) { lr_ = lr; }

    /** @brief 获取学习率 */
    [[nodiscard]] double learning_rate() const { return lr_; }

    /** @brief 获取当前迭代次数 */
    [[nodiscard]] int iteration() const { return iteration_; }

private:
    double lr_;
    int iteration_ = 0;
};

} // namespace vde::sdf