#pragma once
/**
 * @file sdf_gradient.h
 * @brief SDF 梯度计算（数值微分、解析梯度、链式法则）
 *
 * 提供两类梯度计算：
 *
 * 1. **空间梯度** ∇f(p) — SDF 在场中的方向导数（法向量方向）
 *    - 有限差分近似
 *    - 各图元的解析梯度
 *    - CSG 操作的链式法则
 *
 * 2. **参数梯度** ∂f/∂θ — SDF 对形状参数的敏感度
 *    - 有限差分计算 ∂f/∂radius, ∂f/∂height 等
 *    - 用于基于梯度的形状优化
 *
 * ## 梯度约定
 *
 * 对于 SDF，梯度指向距离增加最快的方向。
 * 在表面上 (f=0)，∇f 等于外法向量（单位长度）。
 * 在外部 (f>0)，∇f 指向远离表面的方向。
 * 在内部 (f<0)，∇f 指向最近的表面。
 *
 * @ingroup sdf
 */
#include "vde/sdf/sdf_primitives.h"
#include "vde/sdf/sdf_operations.h"
#include "vde/core/point.h"
#include <cmath>
#include <functional>

namespace vde::sdf {

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

// ═══════════════════════════════════════════════
//  Finite-Difference Gradient
// ═══════════════════════════════════════════════

/**
 * @brief SDF 的空间梯度（中心差分近似）
 *
 * 使用三轴中心有限差分:
 * ∂f/∂x ≈ (f(x+h, y, z) - f(x-h, y, z)) / (2h)
 *
 * 二阶精度，适用于任意黑盒 SDF 函数。
 *
 * @param f SDF 函数（接受 Point3D 参数）
 * @param p 求值点
 * @param h 差分步长（默认 1e-6）
 * @return 梯度向量 ∇f(p)
 *
 * @note 步长太小会导致舍入误差主导，太大导致截断误差主导。
 *       默认值 1e-6 在单精度场景下是较好的平衡。
 * @note 对于含不连续性的 SDF（如 sharp CSG），有限差分可能在边界产生错误梯度。
 *
 * @code{.cpp}
 * auto f = [](const Point3D& p) { return sphere(p, 1.0); };
 * Vector3D n = gradient(f, Point3D(1, 0, 0));  // ≈ (1, 0, 0)
 * @endcode
 *
 * @see evaluate_with_gradient 一次性获取值和梯度
 */
[[nodiscard]] inline Vector3D gradient(
    const std::function<double(const Point3D&)>& f,
    const Point3D& p, double h = 1e-6)
{
    Vector3D g;
    double inv_2h = 0.5 / h;
    g.x() = (f(Point3D(p.x() + h, p.y(), p.z())) - f(Point3D(p.x() - h, p.y(), p.z()))) * inv_2h;
    g.y() = (f(Point3D(p.x(), p.y() + h, p.z())) - f(Point3D(p.x(), p.y() - h, p.z()))) * inv_2h;
    g.z() = (f(Point3D(p.x(), p.y(), p.z() + h)) - f(Point3D(p.x(), p.y(), p.z() - h))) * inv_2h;
    return g;
}

// ═══════════════════════════════════════════════
//  Gradient Result Struct
// ═══════════════════════════════════════════════

/**
 * @brief 前向模梯度结果: SDF 值 + 梯度
 *
 * 在一次计算中返回 SDF 值和空间梯度，避免重复求值。
 * 梯度向量在表面上即为单位外法向量的近似。
 */
struct GradResult {
    double value;     ///< SDF 值 f(p)
    Vector3D grad;    ///< 空间梯度 ∇f(p) = (df/dx, df/dy, df/dz)
};

/**
 * @brief 同时计算 SDF 值和梯度的便捷函数
 *
 * @param f SDF 函数
 * @param p 求值点
 * @param h 差分步长
 * @return value + grad 结构
 *
 * @see gradient 仅计算梯度
 */
[[nodiscard]] inline GradResult evaluate_with_gradient(
    const std::function<double(const Point3D&)>& f,
    const Point3D& p, double h = 1e-6)
{
    return GradResult{f(p), gradient(f, p, h)};
}

// ═══════════════════════════════════════════════
//  Parameter Gradients (Finite Difference)
// ═══════════════════════════════════════════════

/**
 * @brief 球体 SDF 对半径的参数梯度 ∂f/∂r
 *
 * 用中心差分计算: (sphere(p, r+h) - sphere(p, r-h)) / (2h)
 *
 * @param p 采样点
 * @param radius 球体半径
 * @param h 差分步长（默认 1e-6）
 * @return ∂f/∂radius
 *
 * @note 对于球体，解析梯度为 -1（球面上），但此函数提供与实现无关的数值验证。
 *
 * @see sphere_param_grad 封装版本
 */
[[nodiscard]] inline double sphere_radius_gradient(const Point3D& p, double radius, double h = 1e-6)
{
    return (sphere(p, radius + h) - sphere(p, radius - h)) / (2.0 * h);
}

/**
 * @brief 长方体 SDF 对各轴半尺寸的参数梯度 ∂f/∂extent[i]
 *
 * @param p 采样点
 * @param extents 半边长
 * @param axis 轴索引: 0=x, 1=y, 2=z
 * @param h 差分步长
 * @return ∂f/∂extent[axis]
 *
 * @see box_param_grad 封装版本
 */
[[nodiscard]] inline double box_extent_gradient(const Point3D& p, const Point3D& extents,
                                                  int axis, double h = 1e-6)
{
    Point3D ext_plus = extents;
    Point3D ext_minus = extents;
    ext_plus[axis] += h;
    ext_minus[axis] -= h;
    return (box(p, ext_plus) - box(p, ext_minus)) / (2.0 * h);
}

/**
 * @brief 圆环 SDF 对主半径的参数梯度 ∂f/∂major_r
 *
 * @param p 采样点
 * @param major_r 主半径
 * @param minor_r 副半径
 * @param h 差分步长
 * @return 参数梯度
 */
[[nodiscard]] inline double torus_major_gradient(const Point3D& p, double major_r, double minor_r,
                                                   double h = 1e-6)
{
    return (torus(p, major_r + h, minor_r) - torus(p, major_r - h, minor_r)) / (2.0 * h);
}

/**
 * @brief 圆环 SDF 对副半径的参数梯度 ∂f/∂minor_r
 *
 * @return 参数梯度
 */
[[nodiscard]] inline double torus_minor_gradient(const Point3D& p, double major_r, double minor_r,
                                                   double h = 1e-6)
{
    return (torus(p, major_r, minor_r + h) - torus(p, major_r, minor_r - h)) / (2.0 * h);
}

/**
 * @brief 圆柱 SDF 对高度的参数梯度 ∂f/∂height
 *
 * @param p 采样点
 * @param radius 圆柱半径
 * @param height 圆柱高度
 * @param h 差分步长
 * @return 参数梯度
 */
[[nodiscard]] inline double cylinder_height_gradient(const Point3D& p, double radius, double height,
                                                       double h = 1e-6)
{
    return (cylinder(p, radius, height + h) - cylinder(p, radius, height - h)) / (2.0 * h);
}

/**
 * @brief 圆柱 SDF 对半径的参数梯度 ∂f/∂radius
 *
 * @return 参数梯度
 */
[[nodiscard]] inline double cylinder_radius_gradient(const Point3D& p, double radius, double height,
                                                       double h = 1e-6)
{
    return (cylinder(p, radius + h, height) - cylinder(p, radius - h, height)) / (2.0 * h);
}

// ═══════════════════════════════════════════════
//  Analytical Gradients (Exact)
// ═══════════════════════════════════════════════

/**
 * @brief 球体 SDF 的解析梯度
 *
 * **公式:** ∇f = p / |p|
 *
 * 即从球心指向采样点的单位径向向量。
 * 注意：球心处 p = (0,0,0)，梯度退化，返回零向量。
 *
 * @param p 采样点
 * @param radius 球半径（当前未使用，保留以保持接口一致）
 * @return 解析梯度（单位向量）
 *
 * @note 在表面上 (|p| = radius)，梯度即为外法向量。
 */
[[nodiscard]] inline Vector3D sphere_gradient_analytic(const Point3D& p, double /*radius*/)
{
    double n = p.norm();
    if (n < 1e-30) return Vector3D::Zero();
    return p / n;
}

/**
 * @brief 长方体 SDF 的解析梯度
 *
 * **算法:**
 * - 外部: 梯度分量非零当且仅当对应坐标超出 semi-extent，方向向外
 * - 内部: 梯度指向最近的表面（单分量）
 *
 * 注意：长方体 SDF 使用平滑近似，边缘处的梯度在内部/外部是分段常数。
 *
 * @param p 采样点
 * @param half_extents 半边长
 * @return 解析梯度
 */
[[nodiscard]] inline Vector3D box_gradient_analytic(const Point3D& p, const Point3D& half_extents)
{
    Vector3D g = Vector3D::Zero();
    if (std::abs(p.x()) > half_extents.x()) g.x() = (p.x() > 0.0 ? 1.0 : -1.0);
    if (std::abs(p.y()) > half_extents.y()) g.y() = (p.y() > 0.0 ? 1.0 : -1.0);
    if (std::abs(p.z()) > half_extents.z()) g.z() = (p.z() > 0.0 ? 1.0 : -1.0);

    // For interior: gradient points toward nearest face
    if (g.norm() < 1e-30 && half_extents.norm() > 1e-30) {
        // Find the dimension where the point is closest to a face
        Point3D q = Point3D(
            half_extents.x() - std::abs(p.x()),
            half_extents.y() - std::abs(p.y()),
            half_extents.z() - std::abs(p.z())
        );
        // Smallest q component = closest face
        if (q.x() <= q.y() && q.x() <= q.z()) g.x() = (p.x() >= 0.0 ? 1.0 : -1.0);
        else if (q.y() <= q.x() && q.y() <= q.z()) g.y() = (p.y() >= 0.0 ? 1.0 : -1.0);
        else g.z() = (p.z() >= 0.0 ? 1.0 : -1.0);
    }
    return g;
}

/**
 * @brief 圆环 SDF 的解析梯度
 *
 * **SDF:** f = √((√(pₓ²+p𝓏²) - R_major)² + p_y²) - R_minor
 *
 * 梯度通过链式法则推导，将环向距离和 Y 方向分别求导。
 *
 * **特殊情况:** 在 Y 轴上 (pₓ=p𝓏=0)，x 和 z 梯度分量退化为 0。
 *
 * @param p 采样点
 * @param major_r 主半径
 * @param minor_r 副半径
 * @return 解析梯度
 */
[[nodiscard]] inline Vector3D torus_gradient_analytic(const Point3D& p, double major_r, double minor_r)
{
    double rho = std::sqrt(p.x() * p.x() + p.z() * p.z());
    if (rho < 1e-30) {
        // On the Y axis: x and z gradients are 0 by symmetry;
        // y gradient = py / sqrt(major² + py²)
        double inner = std::sqrt(major_r * major_r + p.y() * p.y());
        if (inner < 1e-30) return Vector3D::Zero();
        return Vector3D(0.0, p.y() / inner, 0.0);
    }
    double inner = std::sqrt((rho - major_r) * (rho - major_r) + p.y() * p.y());
    if (inner < 1e-30) return Vector3D::Zero();
    double d_rho = (rho - major_r) / inner;
    return Vector3D(
        d_rho * p.x() / rho,
        p.y() / inner,
        d_rho * p.z() / rho
    );
}

/**
 * @brief 无限圆柱（沿 Y 轴）的解析梯度
 *
 * **公式:** ∇f = (pₓ, 0, p𝓏) / √(pₓ² + p𝓏²)
 *
 * 梯度仅含水平分量，方向为从轴线径向向外。
 * 在轴线上退化返回 (0,1,0)。
 *
 * @param p 采样点
 * @param radius 圆柱半径（未使用，保留接口一致）
 * @return 解析梯度
 */
[[nodiscard]] inline Vector3D cylinder_gradient_analytic(const Point3D& p, double /*radius*/)
{
    double rho = std::sqrt(p.x() * p.x() + p.z() * p.z());
    if (rho < 1e-30) return Vector3D::UnitY(); // degenerate case: along the axis
    return Vector3D(p.x() / rho, 0.0, p.z() / rho);
}

/**
 * @brief 平面的解析梯度
 *
 * **公式:** ∇f = normal
 *
 * 平面的梯度为常向量（法向量），与位置无关。
 *
 * @param normal 单位法向量
 * @return 梯度（等于 normal）
 *
 * @note 这是最简单的解析梯度，因为平面 SDF 是线性的。
 */
[[nodiscard]] inline Vector3D plane_gradient_analytic(const Vector3D& normal)
{
    return normal;
}

/**
 * @brief 胶囊体 SDF 的解析梯度
 *
 * 梯度指向点 p 到最近线段上的投影方向。
 *
 * **算法:**
 * 1. 找到 p 在线段 ab 上的投影点 closest
 * 2. 若 |p-closest| ≈ 0（点在轴上），返回零向量
 * 3. 否则梯度 = (p - closest) / |p - closest|
 *
 * 退化情况: 当 a ≈ b（线段退化为点）时，退化为球体的径向梯度。
 *
 * @param p 采样点
 * @param a 线段起点
 * @param b 线段终点
 * @param radius 扫掠球半径
 * @return 解析梯度
 */
[[nodiscard]] inline Vector3D capsule_gradient_analytic(const Point3D& p,
                                                          const Point3D& a, const Point3D& b,
                                                          double radius)
{
    Point3D pa = p - a;
    Point3D ba = b - a;
    double ba_sq = ba.squaredNorm();
    if (ba_sq < 1e-30) {
        // Degenerate to sphere
        double n = pa.norm();
        if (n < 1e-30) return Vector3D::Zero();
        return pa / n;
    }
    double h = std::clamp(pa.dot(ba) / ba_sq, 0.0, 1.0);
    Point3D closest = a + ba * h;
    Point3D delta = p - closest;
    double d = delta.norm();
    if (d < 1e-30) return Vector3D::Zero();
    return delta / d;
}

// ═══════════════════════════════════════════════
//  Chain Rule for CSG Operations
// ═══════════════════════════════════════════════

/**
 * @brief 并集操作后的梯度传递
 *
 * 取两个子对象中较近者的梯度（因为距离值来自较近者）。
 *
 * @param grad_a 子对象 A 的梯度
 * @param grad_b 子对象 B 的梯度
 * @param d1 A 的 SDF 值
 * @param d2 B 的 SDF 值
 * @return 合并后的梯度
 *
 * @warning 当 d1 ≈ d2 时，梯度可能在两个子对象间不连续跳变。
 *
 * @see chain_smooth_union 避免不连续的平滑版本
 */
[[nodiscard]] inline Vector3D chain_union(const Vector3D& grad_a, const Vector3D& grad_b,
                                            double d1, double d2)
{
    return (d1 < d2) ? grad_a : grad_b;
}

/**
 * @brief 交集操作后的梯度传递
 *
 * 取两个子对象中较远者的梯度（因为距离值来自较远者）。
 *
 * @return 合并后的梯度
 */
[[nodiscard]] inline Vector3D chain_intersection(const Vector3D& grad_a, const Vector3D& grad_b,
                                                   double d1, double d2)
{
    return (d1 > d2) ? grad_a : grad_b;
}

/**
 * @brief 差集操作后的梯度传递
 *
 * 若 A 决定结果 (d1 > -d2)，使用 ∇A；否则使用 -∇B（B 的内部变成外部）。
 *
 * @param grad_a 子对象 A 的梯度
 * @param grad_b 子对象 B 的梯度
 * @param d1 f_A
 * @param d2 f_B
 * @return 合并后的梯度
 */
[[nodiscard]] inline Vector3D chain_difference(const Vector3D& grad_a, const Vector3D& grad_b,
                                                 double d1, double d2)
{
    return (d1 > -d2) ? grad_a : Vector3D(-grad_b.x(), -grad_b.y(), -grad_b.z());
}

/**
 * @brief 平滑并集操作的梯度传递（含完整链式法则）
 *
 * 与普通并集的硬切换不同，本函数使用多项式混合函数的导数来平滑插值两个子对象的梯度。
 *
 * **推导:**
 * f  = d₁(1-h) + d₂h - k·h(1-h)，h = clamp(0.5+0.5(d₂-d₁)/k, 0, 1)
 * ∂f/∂d₁ = 2 - 3h
 * ∂f/∂d₂ = 3h - 1
 *
 * 在混合区域 (0 < h < 1) 内，权重平滑地从 A 过渡到 B。
 *
 * @param grad_a 子对象 A 的梯度
 * @param grad_b 子对象 B 的梯度
 * @param d1 f_A
 * @param d2 f_B
 * @param k 混合参数
 * @return 平滑合并后的梯度
 *
 * @see op_smooth_union 对应的平滑并集 SDF 操作
 * @see chain_union 硬切换版本
 */
[[nodiscard]] inline Vector3D chain_smooth_union(const Vector3D& grad_a, const Vector3D& grad_b,
                                                   double d1, double d2, double k)
{
    if (k <= 0.0) return chain_union(grad_a, grad_b, d1, d2);
    double h = std::clamp(0.5 + 0.5 * (d1 - d2) / k, 0.0, 1.0);
    if (h <= 0.0) return grad_a;
    if (h >= 1.0) return grad_b;

    // Full derivatives w.r.t. d1, d2:
    // ∂f/∂d1 = 2 - 3h,  ∂f/∂d2 = 3h - 1
    double w1 = 2.0 - 3.0 * h;
    double w2 = 3.0 * h - 1.0;
    return grad_a * w1 + grad_b * w2;
}

// ═══════════════════════════════════════════════
//  Chain Rule for Domain Transforms
// ═══════════════════════════════════════════════

/**
 * @brief 平移变换后的梯度传递（梯度不变）
 *
 * 平移是等距变换（isometry），函数值不变，
 * 因此子对象的梯度直接穿透。
 *
 * @param child_result 变换空间中的子对象求值结果
 * @param offset 平移量（仅接口保留）
 * @return 相同梯度
 */
[[nodiscard]] inline GradResult chain_translate(const GradResult& child_result, const Point3D& /*offset*/)
{
    return child_result;
}

/**
 * @brief 绕 Y 轴旋转后的梯度传递
 *
 * 由于 op_rotate 对采样点施加 R(-θ)，链式法则要求对梯度施加 R(θ)。
 * 即：∇f_world = R(θ) · ∇f_local
 *
 * @param child_result 旋转后坐标系中的子对象求值
 * @param angle_rad 旋转角度
 * @return 世界坐标系的梯度
 *
 * @see op_rotate 对应的域变形操作
 */
[[nodiscard]] inline GradResult chain_rotate(const GradResult& child_result, double angle_rad)
{
    GradResult result = child_result;
    double c = std::cos(angle_rad);
    double s = std::sin(angle_rad);
    // R(θ) = [[c, 0, -s], [0, 1, 0], [s, 0, c]]
    result.grad = Vector3D(
        c * child_result.grad.x() - s * child_result.grad.z(),
        child_result.grad.y(),
        s * child_result.grad.x() + c * child_result.grad.z()
    );
    return result;
}

/**
 * @brief 非均匀缩放后的梯度传递
 *
 * **公式:** ∇f_world = (∇f_child.x / sx, ∇f_child.y / sy, ∇f_child.z / sz)
 *
 * 因为 f_world(p) = f_child(p / s)，由链式法则得 ∂f_world/∂p = ∇f_child · diag(1/s)。
 *
 * @note 缩放后梯度不再是单位长度，除非是均匀缩放。
 *
 * @param child_result 缩放后坐标系中的子对象求值
 * @param factors 缩放因子
 * @return 世界坐标系的梯度
 */
[[nodiscard]] inline GradResult chain_scale(const GradResult& child_result, const Point3D& factors)
{
    GradResult result = child_result;
    result.grad = Vector3D(
        child_result.grad.x() / factors.x(),
        child_result.grad.y() / factors.y(),
        child_result.grad.z() / factors.z()
    );
    return result;
}

/**
 * @brief 扭曲变换后的梯度传递（含完整雅可比）
 *
 * 扭曲 T(p) 绕 Y 轴旋转角度 amount·p_y。
 * 梯度传递需要对雅可比矩阵的转置做链式法则：
 *
 * **雅可比 J:**
 * ```
 * ∂T_x/∂p_x = c,  ∂T_x/∂p_y = -a·(p_x·s + p_z·c),  ∂T_x/∂p_z = -s
 * ∂T_y/∂p_x = 0,  ∂T_y/∂p_y = 1,                     ∂T_y/∂p_z = 0
 * ∂T_z/∂p_x = s,  ∂T_z/∂p_y =  a·(p_x·c - p_z·s),   ∂T_z/∂p_z = c
 * ```
 * 其中 c = cos(a·p_y), s = sin(a·p_y), a = amount。
 *
 * @param child_result 扭曲空间中的子对象求值
 * @param amount 单位高度的旋转量
 * @param p 原始采样点
 * @return 世界坐标系的梯度
 *
 * @see op_twist 对应的域变形操作
 * @see chain_rotate 不含 Y 依赖的简化旋转梯度传递
 */
[[nodiscard]] inline GradResult chain_twist(const GradResult& child_result,
                                              double amount, const Point3D& p)
{
    double angle = amount * p.y();
    double c = std::cos(angle);
    double s = std::sin(angle);

    // Forward-twisted position
    double tp_x = p.x() * c - p.z() * s;
    double tp_z = p.x() * s + p.z() * c;

    GradResult result = child_result;

    // Spatial part: inverse rotation R(-angle) applied to child gradient
    double gx = c * child_result.grad.x() + s * child_result.grad.z();
    double gz = -s * child_result.grad.x() + c * child_result.grad.z();

    // Y-component: extra term from angle's dependence on p.y
    // ∂/∂y correction: amount * (tp_x * gz_child - tp_z * gx_child)
    double gy = child_result.grad.y() + amount * (tp_x * child_result.grad.z()
                                                 - tp_z * child_result.grad.x());

    result.grad = Vector3D(gx, gy, gz);
    return result;
}

/**
 * @brief 周期重复后的梯度传递（梯度不变）
 *
 * 重复映射是平移变换的叠加，每个晶格内梯度相同（周期域）。
 * 在晶格边界处存在不连续性。
 *
 * @param child_result 中心晶格内的子对象求值
 * @param cell 晶格尺寸（仅接口保留）
 * @param p 原始采样点（仅接口保留）
 * @return 相同梯度
 *
 * @warning 在晶格边界处梯度不连续，可能影响基于梯度的优化。
 */
[[nodiscard]] inline GradResult chain_repeat(const GradResult& child_result,
                                               const Point3D& /*cell*/, const Point3D& /*p*/)
{
    return child_result;
}

// ═══════════════════════════════════════════════
//  Parameter Gradient Vector (All Shape Params)
// ═══════════════════════════════════════════════

/**
 * @brief 参数梯度结构体
 *
 * 存储 SDF 对所有形状参数的偏导数：
 * ∂f/∂radius, ∂f/∂height, ∂f/∂extent[3], ∂f/∂major, ∂f/∂minor,
 * ∂f/∂angle, ∂f/∂translate[3], ∂f/∂rotate, ∂f/∂scale[3]
 *
 * 用于基于梯度的形状优化和参数拟合。
 *
 * @see sphere_param_grad, box_param_grad, cylinder_param_grad
 */
struct ParamGrad {
    double d_radius = 0.0;             ///< ∂f/∂radius
    double d_height = 0.0;             ///< ∂f/∂height
    Point3D d_extents = Point3D::Zero(); ///< ∂f/∂(extent.x, extent.y, extent.z)
    double d_major = 0.0;              ///< ∂f/∂major_radius
    double d_minor = 0.0;              ///< ∂f/∂minor_radius
    double d_angle = 0.0;              ///< ∂f/∂angle_rad
    Vector3D d_normal = Vector3D::Zero(); ///< ∂f/∂normal (无效分量清零)
    double d_offset = 0.0;             ///< ∂f/∂offset
    Point3D d_translate = Point3D::Zero(); ///< ∂f/∂translate_offset
    double d_rotate = 0.0;             ///< ∂f/∂rotation_angle
    Point3D d_scale = Point3D::Zero(); ///< ∂f/∂scale_factors
};

/**
 * @brief 球体所有参数的梯度（当前仅 radius 有意义）
 *
 * @param p 采样点
 * @param radius 球体半径
 * @param h 差分步长
 * @return 参数梯度（仅 d_radius 非零）
 *
 * @see sphere_radius_gradient 底层数值差分
 */
[[nodiscard]] inline ParamGrad sphere_param_grad(const Point3D& p, double radius, double h = 1e-6)
{
    ParamGrad pg;
    pg.d_radius = sphere_radius_gradient(p, radius, h);
    return pg;
}

/**
 * @brief 长方体所有参数的梯度
 *
 * 分别计算三个轴的半尺寸对 SDF 值的影响。
 *
 * @param p 采样点
 * @param half_extents 半边长
 * @param h 差分步长
 * @return 参数梯度（d_extents 分量填充）
 *
 * @see box_extent_gradient 底层数值差分
 */
[[nodiscard]] inline ParamGrad box_param_grad(const Point3D& p, const Point3D& half_extents,
                                                double h = 1e-6)
{
    ParamGrad pg;
    pg.d_extents.x() = box_extent_gradient(p, half_extents, 0, h);
    pg.d_extents.y() = box_extent_gradient(p, half_extents, 1, h);
    pg.d_extents.z() = box_extent_gradient(p, half_extents, 2, h);
    return pg;
}

/**
 * @brief 圆柱所有参数的梯度
 *
 * 分别计算半径和高度对 SDF 值的影响。
 *
 * @param p 采样点
 * @param radius 圆柱半径
 * @param height 圆柱高度
 * @param h 差分步长
 * @return 参数梯度（d_radius, d_height 填充）
 */
[[nodiscard]] inline ParamGrad cylinder_param_grad(const Point3D& p, double radius, double height,
                                                     double h = 1e-6)
{
    ParamGrad pg;
    pg.d_radius = cylinder_radius_gradient(p, radius, height, h);
    pg.d_height = cylinder_height_gradient(p, radius, height, h);
    return pg;
}

} // namespace vde::sdf