#pragma once
/**
 * @file sdf_primitives.h
 * @brief 有符号距离函数（SDF）基础图元
 *
 * 提供常见几何体到三维空间中各点的有符号距离计算。
 * 所有图元定义在 vde::sdf 命名空间下，使用 vde::core::Point3D 作为输入。
 *
 * ## SDF 约定
 *
 * - f(p) < 0: 点 p 在几何体**内部**
 * - f(p) = 0: 点 p 在几何体**表面**上
 * - f(p) > 0: 点 p 在几何体**外部**
 * - |f(p)| 的绝对值近似为 p 到最近表面的欧几里得距离
 *
 * ## 数值特性
 *
 * - 所有内联函数标记为 `[[nodiscard]]`，强制调用方使用返回值
 * - 距离函数在表面附近满足 |∇f| ≈ 1（Lipschitz 连续，步进可靠）
 * - 部分复杂图元（cone, triangular_prism, link, infinite_cone）为非内联实现，定义在 src/sdf/ 中
 *
 * @note 所有图元假设未施加变换——坐标系原点即图元中心/基点。
 *       需要平移/旋转/缩放时，使用 sdf_operations.h 中的域变形算子对采样点预先变换。
 *
 * @ingroup sdf
 */
#include "vde/core/point.h"
#include <cmath>
#include <algorithm>

namespace vde::sdf {
using core::Point3D;
using core::Vector3D;

// ═══════════════════════════════════════════════
//  Basic Primitives
// ═══════════════════════════════════════════════

/**
 * @brief 球体的 SDF
 *
 * 计算三维空间中一点到球面的有符号距离。
 *
 * **数学公式:** f(p) = |p| - r
 *
 * ```
 * f < 0: 点在球内
 * f = 0: 点在球面上
 * f > 0: 点在球外
 * ```
 *
 * |∇f| = 1 在球外和球面上严格成立，球内在原点处退化（梯度为零）。
 *
 * @param p 采样点（三维坐标）
 * @param radius 球半径，必须 > 0
 * @return 有符号距离
 *
 * @note 球位于原点。如需其他位置，使用 op_translate 对采样点做反向平移。
 *
 * @code{.cpp}
 * double d = sphere(Point3D(1.0, 0.0, 0.0), 0.5);  // d = 0.5 (外部)
 * double inside = sphere(Point3D(0.0, 0.0, 0.0), 1.0); // d = -1.0 (球心)
 * @endcode
 *
 * @see ellipsoid 各轴半径不同的椭球
 * @see cylinder 沿 Y 轴的圆柱
 */
[[nodiscard]] inline double sphere(const Point3D& p, double radius) {
    return p.norm() - radius;
}

/**
 * @brief 轴对齐立方体的 SDF
 *
 * 计算点到以原点为中心的轴对齐长方体表面的有符号距离。
 *
 * **算法概要:**
 * 1. 将点的各分量取绝对值，得到第一卦限对称映射
 * 2. 减去半边长得到相对于面的距离 qᵢ = |pᵢ| - halfExtentᵢ
 * 3. 外部距离 = ||max(q,0)||，内部惩罚 = min(max(qₓ,q_y,q_z), 0)
 *
 * 结果在面外为真实欧氏距离，边/角区域使用 L2 距离组合保证连续性。
 *
 * @param p 采样点
 * @param half_extents 半边长 (hx, hy, hz)，均必须 ≥ 0
 * @return 有符号距离
 *
 * @code{.cpp}
 * // 2×3×4 的盒子
 * double d = box(Point3D(0.0, 0.0, 0.0), Point3D(1.0, 1.5, 2.0));
 * @endcode
 *
 * @see round_box 带圆角的盒子
 * @see wedge 对角切割的半盒子
 */
[[nodiscard]] inline double box(const Point3D& p, const Point3D& half_extents) {
    Point3D q = p.cwiseAbs() - half_extents;
    return q.cwiseMax(0.0).norm()
         + std::min(std::max({q.x(), q.y(), q.z()}), 0.0);
}

/**
 * @brief 圆角立方体的 SDF
 *
 * 与 box() 使用相同的核心算法，但在最后减去圆角半径 r。
 * 相当于先对 box SDF 做膨胀（offset），再缩回，等价于对边角做半径为 r 的倒圆处理。
 *
 * **数学公式:** f(p) = box(p, half_extents) - r
 *
 * @param p 采样点
 * @param half_extents 半边长 (hx, hy, hz)
 * @param r 圆角半径，必须 < min(half_extents) 否则整体形态改变
 * @return 有符号距离
 *
 * @warning 当 r 超过某轴半边长时，该维度的边会被完全圆化，表面不再平坦。
 *          实际应用中限制 r ≤ 0.5 * min(hx, hy, hz) 可保证形态可控。
 *
 * @code{.cpp}
 * // 带 0.2 圆角的 2×2×2 立方体
 * double d = round_box(p, Point3D(1, 1, 1), 0.2);
 * @endcode
 *
 * @see box 无圆角的长方体
 * @see op_round 对任意 SDF 做圆角偏移
 */
[[nodiscard]] inline double round_box(const Point3D& p, const Point3D& half_extents, double r) {
    Point3D q = p.cwiseAbs() - half_extents;
    return q.cwiseMax(0.0).norm()
         + std::min(std::max({q.x(), q.y(), q.z()}), 0.0)
         - r;
}

/**
 * @brief 圆环面的 SDF
 *
 * 圆环位于 XZ 平面内，绕 Y 轴对称。
 *
 * **数学公式:**
 * 1. 计算环向半径: q = √(pₓ² + p𝓏²) - R_major
 * 2. 管向距离: f = √(q² + p_y²) - R_minor
 *
 * R_major = 环心到管心的距离（主半径）
 * R_minor = 管的截面半径（副半径）
 *
 * @param p 采样点
 * @param major_radius 主半径（环的半径），必须 > 0
 * @param minor_radius 副半径（管的粗细），必须 > 0
 * @return 有符号距离
 *
 * @warning 当 major_radius ≤ minor_radius 时，环孔消失形成苹果状（self-intersecting torus）。
 *
 * @code{.cpp}
 * // 主半径 2.0、管半径 0.5 的圆环
 * double d = torus(Point3D(2.0, 0.0, 0.0), 2.0, 0.5);  // d ≈ 0.0（表面上）
 * @endcode
 *
 * @see link 双环连接体（两个并行圆环的组合）
 */
[[nodiscard]] inline double torus(const Point3D& p, double major_radius, double minor_radius) {
    double qx = std::sqrt(p.x() * p.x() + p.z() * p.z()) - major_radius;
    return std::sqrt(qx * qx + p.y() * p.y()) - minor_radius;
}

/**
 * @brief 胶囊体的 SDF
 *
 * 以线段 ab 为骨架、半径为 radius 的扫掠球体。
 * 等价于线段上最近的点到 p 的距离减去 radius。
 *
 * **算法概要:**
 * 1. 将 p 投影到线段 ab 上，投影参数 h = clamp(pa·ba / ba·ba, 0, 1)
 * 2. 计算 p 到投影点的距离减去 radius
 *
 * 退化情况: 当 |ba| ≈ 0（a ≈ b），退化为以 a 为球心的球体。
 *
 * @param p 采样点
 * @param a 线段起点
 * @param b 线段终点
 * @param radius 扫掠球半径，必须 > 0
 * @return 有符号距离
 *
 * @code{.cpp}
 * // 从 (0,-2,0) 到 (0,2,0) 的胶囊、半径 0.5
 * double d = capsule(p, Point3D(0,-2,0), Point3D(0,2,0), 0.5);
 * @endcode
 *
 * @see cylinder 两端平头的圆柱（非球形端盖）
 */
[[nodiscard]] inline double capsule(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
        return pa.norm() - radius;
    }
    double h = std::clamp(pa.dot(ba) / ba_sq, 0.0, 1.0);
    return (pa - ba * h).norm() - radius;
}

/**
 * @brief 带平端盖的圆柱体 SDF
 *
 * 沿 Y 轴、以原点为中心的有限长圆柱体。两端为平切面（非球形端盖）。
 *
 * **算法概要:**
 * 1. 水平距离: d_xz = √(pₓ² + p𝓏²) - radius
 * 2. 垂直距离: d_y = |p_y| - height/2
 * 3. 用二维 min/max L2 组合计算矩形区域到圆角的 SDF
 *
 * 结果在柱面上为真实距离，在端盖边缘有连续过渡。
 *
 * @param p 采样点
 * @param radius 圆柱截面半径，必须 > 0
 * @param height 圆柱总高度（从 -h/2 到 +h/2），必须 > 0
 * @return 有符号距离
 *
 * @code{.cpp}
 * // 半径 1.0、高度 3.0 的圆柱
 * double d = cylinder(Point3D(1.0, 0.0, 0.0), 1.0, 3.0);
 * @endcode
 *
 * @see capsule 两端球形端盖的圆柱
 * @see infinite_cylinder 无限长圆柱
 * @see cone 圆锥台
 */
[[nodiscard]] inline double cylinder(const Point3D& p, double radius, double height) {
    double d_xz = std::sqrt(p.x() * p.x() + p.z() * p.z()) - radius;
    double d_y  = std::abs(p.y()) - height * 0.5;
    return std::min(std::max(d_xz, d_y), 0.0)
         + std::sqrt(std::max(d_xz, 0.0) * std::max(d_xz, 0.0)
                   + std::max(d_y,  0.0) * std::max(d_y,  0.0));
}

/**
 * @brief 平面的 SDF
 *
 * 通过原点的平面，法向量为 normal，偏移量为 offset。
 *
 * **数学公式:** f(p) = p·n - offset
 *
 * 其中 n 为单位法向量。f > 0 表示点在法向量指向的一侧。
 * 偏移量沿法向量正方向移动平面。
 *
 * @param p 采样点
 * @param normal 单位法向量（必须已归一化）
 * @param offset 沿法向量正方向的偏移量
 * @return 有符号距离
 *
 * @pre normal.norm() ≈ 1.0（调用方负责归一化）
 *
 * @code{.cpp}
 * // 通过原点、法向量朝上的平面
 * double d = plane(p, Vector3D(0,1,0), 0.0);
 * // 向上偏移 2 个单位的平面
 * double d2 = plane(p, Vector3D(0,1,0), 2.0);
 * @endcode
 *
 * @see wedge 两个平面+盒子的组合（对角切割）
 */
[[nodiscard]] inline double plane(const Point3D& p, const Vector3D& normal, double offset) {
    return p.dot(normal) - offset;
}

/**
 * @brief 椭球体的 SDF（有界梯度近似）
 *
 * 以原点为中心的椭球，半轴长为 radii = (rx, ry, rz)。
 *
 * **算法概要:**
 *
 * 使用 bound-normalized 近似方法：
 * 1. 缩放坐标: pr = p / radii
 * 2. k0 = |pr|，k1 = |p / radii²|
 * 3. f = k0 * (k0 - 1) / k1
 *
 * 该方法保证在表面上 f = 0，但内部/远处的距离为近似值。
 * 相比精确椭球距离，避开了六次方程的昂贵求根。
 *
 * @param p 采样点
 * @param radii 半轴长 (rx, ry, rz)，均必须 > 0
 * @return 有符号距离（近似）
 *
 * @note 当三个半径相等时退化为球体，此时使用 sphere() 更精确高效。
 * @note 当点到原点距离远大于最大半径时，近似精度下降。
 *
 * @code{.cpp}
 * // 扁椭球: rx=2, ry=1, rz=1
 * double d = ellipsoid(Point3D(2.0, 0.0, 0.0), Point3D(2, 1, 1));
 * @endcode
 *
 * @see sphere 精确球体（rx=ry=rz 时的特例）
 */
[[nodiscard]] inline double ellipsoid(const Point3D& p, const Point3D& radii) {
    Point3D pr(p.x() / radii.x(), p.y() / radii.y(), p.z() / radii.z());
    double k0 = pr.norm();
    double k1 = Point3D(p.x() / (radii.x() * radii.x()),
                         p.y() / (radii.y() * radii.y()),
                         p.z() / (radii.z() * radii.z())).norm();
    if (k1 < 1e-30) {
        // Degenerate: point at origin → distance = -min(radii)
        return -std::min({radii.x(), radii.y(), radii.z()});
    }
    return k0 * (k0 - 1.0) / k1;
}

// ═══════════════════════════════════════════════
//  Parametric Primitives
// ═══════════════════════════════════════════════

/**
 * @brief 正六棱柱的 SDF
 *
 * 六边形位于 XZ 平面，尖顶沿 X 轴方向，沿 Y 轴挤出的正六棱柱。
 *
 * **算法概要:**
 * 1. 将点对称映射到第一象限后，计算二维六边形的 SDF
 * 2. 六边形 SDF 通过两组约束平面表示: max(px + 0.577*pz, 1.155*pz) - radius
 * 3. 沿 Y 轴做 bounded extrusion
 *
 * @param p 采样点
 * @param radius 六边形外接圆半径（中心到顶点），必须 > 0
 * @param height 挤出总高度（±h/2 沿 Y 轴），必须 > 0
 * @return 有符号距离
 *
 * @code{.cpp}
 * // 外接圆半径 1.0、高度 2.0 的六棱柱
 * double d = hex_prism(Point3D(0.0, 0.0, 0.0), 1.0, 2.0);
 * @endcode
 *
 * @see extrusion_bounded 通用 bounded extrusion 算子
 */
[[nodiscard]] inline double hex_prism(const Point3D& p, double radius, double height) {
    double px = std::abs(p.x());
    double pz = std::abs(p.z());
    // 2D hexagon SDF (pointy-top in XZ plane)
    double d_hex = std::max(px + pz * 0.577350269189626, pz * 1.154700538379252) - radius;
    // Extrude along Y
    return std::max(d_hex, std::abs(p.y()) - height * 0.5);
}

/**
 * @brief 无限长圆柱的 SDF
 *
 * 沿任意方向的无限长圆柱（无端盖），轴线过原点。
 *
 * **数学公式:** f(p) = |p × axis| - radius
 *
 * 叉积给出点到轴线的垂直距离。
 *
 * @param p 采样点
 * @param axis 圆柱轴线方向（必须归一化，|axis| = 1）
 * @param radius 截面半径，必须 > 0
 * @return 有符号距离
 *
 * @pre axis.norm() ≈ 1.0
 *
 * @code{.cpp}
 * // 沿 Z 轴的无限圆柱
 * double d = infinite_cylinder(p, Vector3D(0,0,1), 0.5);
 * @endcode
 *
 * @see cylinder 有限长圆柱（带平端盖）
 */
[[nodiscard]] inline double infinite_cylinder(const Point3D& p, const Vector3D& axis, double radius) {
    return p.cross(axis).norm() - radius;
}

/**
 * @brief 楔形体的 SDF
 *
 * 立方体沿 XZ 平面对角线切割的一半。切割面为 z = x。
 * 楔形体占据盒子内 z ≥ x 的区域。
 *
 * **几何形状:**
 * - 盒子: [-w/2, w/2] × [-h/2, h/2] × [-d/2, d/2]
 * - 保留区域: 盒子内 z ≥ x 的半空间
 *
 * @param p 采样点
 * @param w 宽度（X 方向）
 * @param h 高度（Y 方向）
 * @param d 深度（Z 方向）
 * @return 有符号距离
 *
 * @code{.cpp}
 * // 2×2×2 的楔形体（保留 z ≥ x 的一半）
 * double d = wedge(Point3D(-0.5, 0.0, 0.5), 2.0, 2.0, 2.0);
 * @endcode
 *
 * @see box 完整立方体
 * @see plane 独立平面 SDF
 */
[[nodiscard]] inline double wedge(const Point3D& p, double w, double h, double d) {
    Point3D q = p.cwiseAbs() - Point3D(w * 0.5, h * 0.5, d * 0.5);
    double d_box = q.cwiseMax(0.0).norm()
                 + std::min(std::max({q.x(), q.y(), q.z()}), 0.0);
    // Diagonal plane x - z ≤ 0  (i.e., z ≥ x)
    double d_cut = (p.x() - p.z()) * 0.707106781186548;  // 1/√2
    return std::max(d_box, d_cut);
}

// ═══════════════════════════════════════════════
//  2D Extrusions & Revolution
// ═══════════════════════════════════════════════

/**
 * @brief 二维 SDF 沿 Z 轴的无限挤出
 *
 * 将二维 SDF 沿 Z 轴无限延伸。三维点的 SDF 值等于其 (x,y) 投影处的二维 SDF 值。
 *
 * **使用模式:** 先计算二维 SDF 值，再传入此函数:
 *
 * @code{.cpp}
 * double sdf_2d = circle_sdf(p.x(), p.y());  // 用户自定义二维 SDF
 * double sdf_3d = extrusion(p, sdf_2d);
 * @endcode
 *
 * @param p 三维采样点（仅 x,y 分量被隐式使用——通过 sdf_2d 参数间接）
 * @param sdf_2d 预计算的二维 SDF 值（在 (p.x, p.y) 处求值）
 * @return 三维有符号距离（等于 sdf_2d）
 *
 * @note p 参数仅保持接口一致，实际不影响结果。
 * @note 对于有限长度的挤出，使用 extrusion_bounded()。
 *
 * @see extrusion_bounded 有限长度的挤出
 * @see revolution 二维轮廓绕 Y 轴旋转
 */
[[nodiscard]] inline double extrusion(const Point3D& /*p*/, double sdf_2d) {
    return sdf_2d;
}

/**
 * @brief 二维 SDF 的有限长度挤出
 *
 * 将二维 SDF 沿 Z 轴挤出并通过 |pz| ≤ half_height 做 bounded 限制。
 * 相当于对二维 SDF 与 Z 轴平板做交集。
 *
 * @param p 三维采样点
 * @param sdf_2d 预计算的二维 SDF 值（在 (p.x, p.y) 处求值）
 * @param half_height 挤出半高度（Z 方向），必须 ≥ 0
 * @return 有符号距离
 *
 * @code{.cpp}
 * // 圆形截面、总高 4.0 的短柱
 * double d2d = std::sqrt(p.x()*p.x() + p.y()*p.y()) - 1.0;
 * double d3d = extrusion_bounded(p, d2d, 2.0);
 * @endcode
 *
 * @see extrusion 无限挤出
 * @see hex_prism 六棱柱（特例化 bounded extrusion 实现）
 */
[[nodiscard]] inline double extrusion_bounded(const Point3D& p, double sdf_2d, double half_height) {
    return std::max(sdf_2d, std::abs(p.z()) - half_height);
}

/**
 * @brief 二维轮廓绕 Y 轴的旋转体 SDF
 *
 * 将二维 SDF 绕 Y 轴旋转，生成旋转体的三维 SDF。
 *
 * **使用模式:**
 * 二维 SDF 定义在 (r, y) 坐标中，其中 r 为水平距离 = √(pₓ² + p𝓏²) - offset。
 * 用户需要在调用前将 (r, y) 转换为二维采样并计算 SDF 值。
 *
 * @param p 三维采样点（接口保留，不直接使用）
 * @param sdf_2d 预计算的二维 SDF（在 (√(pₓ²+p𝓏²)-offset, p.y) 处求值）
 * @param offset 径向偏移量（从 Y 轴的起点距离）
 * @return 有符号距离
 *
 * @code{.cpp}
 * // 在距 Y 轴 2.0 处、截面半径 0.5 的圆环
 * double r = std::sqrt(p.x()*p.x() + p.z()*p.z()) - 2.0;
 * double sdf_2d = std::sqrt(r*r + p.y()*p.y()) - 0.5;
 * double d3d = revolution(p, sdf_2d, 2.0);
 * @endcode
 *
 * @see torus 圆环（revolution 的特例：截面为圆）
 */
[[nodiscard]] inline double revolution(const Point3D& /*p*/, double sdf_2d, double /*offset*/) {
    return sdf_2d;
}

// ═══════════════════════════════════════════════
//  Complex Primitives (implemented in src/sdf/)
// ═══════════════════════════════════════════════

/**
 * @brief 圆锥台的 SDF（非内联实现）
 *
 * 沿 Y 轴、以原点为中心的有限长圆锥台。
 * 顶点在 y = +height/2（半径 0），底面在 y = -height/2（半径 height·tan(angle_rad)）。
 *
 * **算法概要:**
 * 将三维问题分解为二维 (r, y) 平面内点到三角形区域的距离计算，
 * 其中 r = √(pₓ² + p𝓏²)，三角形由底边和高决定。
 *
 * @param p 采样点
 * @param angle_rad 半锥角（弧度），决定底面半径
 * @param height 锥体总高度，必须 > 0
 * @return 有符号距离
 *
 * @warning 当 angle_rad 接近 π/2 时底面半径极大，可能导致数值不稳定。
 *
 * @see cylinder 半锥角为 0 时的退化情况
 * @see infinite_cone 无限圆锥
 */
[[nodiscard]] double cone(const Point3D& p, double angle_rad, double height);

/**
 * @brief 三角棱柱的 SDF（非内联实现）
 *
 * 以 XY 平面中任意三角形为截面、沿 Z 轴 ±height/2 挤出的棱柱。
 *
 * @param p 采样点
 * @param a 三角形顶点 A（XY 平面内）
 * @param b 三角形顶点 B（XY 平面内）
 * @param c 三角形顶点 C（XY 平面内）
 * @param height 挤出总高度（沿 Z 轴），必须 > 0
 * @return 有符号距离
 *
 * @pre a, b, c 三点不共线（构成有效三角形）
 *
 * @code{.cpp}
 * // 顶点为 (0,0), (2,0), (0,2)、高度 1.0 的三棱柱
 * double d = triangular_prism(p, Point3D(0,0,0), Point3D(2,0,0), Point3D(0,2,0), 1.0);
 * @endcode
 *
 * @see wedge 直角三角形楔形体
 */
[[nodiscard]] double triangular_prism(const Point3D& p, const Point3D& a, const Point3D& b,
                                      const Point3D& c, double height);

/**
 * @brief 双环链接体的 SDF（非内联实现）
 *
 * 两个并列的圆环面沿 X 轴排列，等效于经典链环结构。
 * 常用于链条、连杆等机械连接件的快速建模。
 *
 * @param p 采样点
 * @param length 两个环心沿 X 轴的距离
 * @param major_r 每个环的主半径（环的半径）
 * @param minor_r 每个环的管半径（截面的粗细）
 * @return 有符号距离
 *
 * @code{.cpp}
 * // 环间距 2.0、主半径 1.0、管半径 0.2 的双环
 * double d = link(Point3D(1.0, 0.0, 0.0), 2.0, 1.0, 0.2);
 * @endcode
 *
 * @see torus 单个圆环
 */
[[nodiscard]] double link(const Point3D& p, double length, double major_r, double minor_r);

/**
 * @brief 无限圆锥的 SDF（非内联实现）
 *
 * 从顶点出发沿轴方向无限延伸的圆锥，半锥角为 angle_rad。
 * 与 cone() 不同，本函数无高度限制，锥体向轴的正方向无限扩展。
 *
 * @param p 采样点
 * @param apex 锥顶位置
 * @param axis 锥轴方向（必须归一化），锥体沿此方向扩展
 * @param angle_rad 半锥角（弧度）
 * @return 有符号距离
 *
 * @pre axis.norm() ≈ 1.0
 *
 * @code{.cpp}
 * // 顶点在原点、沿 Y 轴、半锥角 30° 的无限锥
 * double d = infinite_cone(p, Point3D(0,0,0), Vector3D(0,1,0), M_PI/6);
 * @endcode
 *
 * @see cone 有限高圆锥台
 */
[[nodiscard]] double infinite_cone(const Point3D& p, const Point3D& apex,
                                   const Vector3D& axis, double angle_rad);

} // namespace vde::sdf