#pragma once
/**
 * @file sdf_operations.h
 * @brief SDF 布尔运算、修饰算子与域变形
 *
 * 本文件提供三类操作：
 *
 * 1. **布尔运算（CSG）** — 在距离值上执行组合逻辑
 * 2. **修饰算子（Modifiers）** — 对距离值做数值变换
 * 3. **域变形（Domain Deformation）** — 对采样点空间做几何变换
 *
 * 所有操作遵循 SDF 惯例：对采样点或其距离值做变换后，保持 f<0=内部、f>0=外部的语义。
 * 域变形算子应用到采样点（而非对象），因此位移/旋转方向与直觉相反（逆变换）。
 *
 * @ingroup sdf
 */
#include "vde/core/point.h"
#include <cmath>
#include <algorithm>

namespace vde::sdf {

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

// ── Boolean operations on distance values ──

/**
 * @brief 布尔并集: A ∪ B
 *
 * 取两个距离值的最小值。在两个子对象的距离场中选择较近者。
 *
 * **数学公式:** f = min(d1, d2)
 *
 * 最常见的 CSG 组合操作，用于合并两个几何体。
 *
 * @param d1 第一个对象的 SDF 值
 * @param d2 第二个对象的 SDF 值
 * @return 并集的 SDF 值
 *
 * @code{.cpp}
 * double d = op_union(sphere(p, 1.0), box(p, Point3D(0.5, 0.5, 0.5)));
 * @endcode
 *
 * @see op_smooth_union 带平滑过渡的并集
 * @see op_intersection 交集
 */
[[nodiscard]] inline double op_union(double d1, double d2) {
    return std::min(d1, d2);
}

/**
 * @brief 布尔交集: A ∩ B
 *
 * 取两个距离值的最大值。保留同时属于两个对象的区域。
 *
 * **数学公式:** f = max(d1, d2)
 *
 * @param d1 第一个对象的 SDF 值
 * @param d2 第二个对象的 SDF 值
 * @return 交集的 SDF 值
 *
 * @code{.cpp}
 * // 球体与盒子的交集（取共同区域）
 * double d = op_intersection(sphere(p, 1.0), box(p, Point3D(0.5, 0.5, 0.5)));
 * @endcode
 *
 * @see op_smooth_intersection 带平滑过渡的交集
 */
[[nodiscard]] inline double op_intersection(double d1, double d2) {
    return std::max(d1, d2);
}

/**
 * @brief 布尔差集: A \ B
 *
 * 从对象 A 中减去对象 B 所占区域。
 *
 * **数学公式:** f = max(d1, -d2)
 *
 * 将第二个距离值取反后取最大值：B 的外部变成内部，内部变成外部。
 *
 * @param d1 被减对象 A 的 SDF 值
 * @param d2 减去对象 B 的 SDF 值
 * @return 差集的 SDF 值
 *
 * @code{.cpp}
 * // 从球体中减去盒子
 * double d = op_difference(sphere(p, 1.0), box(p, Point3D(0.3, 0.3, 0.3)));
 * @endcode
 *
 * @see op_smooth_difference 带平滑过渡的差集
 */
[[nodiscard]] inline double op_difference(double d1, double d2) {
    return std::max(d1, -d2);
}

/**
 * @brief 平滑并集（带混合过渡）
 *
 * 通过多项式混合函数在两个对象接合处产生平滑圆角过渡。
 *
 * **数学公式:**
 * h = clamp(0.5 + 0.5·(d1-d2)/k, 0, 1)
 * f = d1·(1-h) + d2·h - k·h·(1-h)
 *
 * 混合参数 k 控制过渡区域的宽度。k 越大，过渡越平滑但越偏离原始几何。
 *
 * @param d1 第一个对象的 SDF 值
 * @param d2 第二个对象的 SDF 值
 * @param k 混合强度参数（> 0 时平滑过渡；≤ 0 时退化为普通并集）
 * @return 平滑混合后的 SDF 值
 *
 * @note 多项式混合是函数 f(x)=x 在 [0,k] 区间的平滑近似。在 d1=d2 处，f 向下偏移 k/4。
 * @warning k 值过大（超过对象尺寸的 10%）会严重改变几何形态。
 *
 * @code{.cpp}
 * // 两个球体平滑融合，过渡区宽度 0.3
 * double d = op_smooth_union(sphere(p, 1.0), sphere(p - Point3D(1.2,0,0), 1.0), 0.3);
 * @endcode
 *
 * @see op_union 无平滑的标准并集
 * @see op_smooth_intersection 平滑交集
 * @see op_smooth_difference 平滑差集
 */
[[nodiscard]] inline double op_smooth_union(double d1, double d2, double k) {
    if (k <= 0.0) return std::min(d1, d2);
    double h = std::clamp(0.5 + 0.5 * (d1 - d2) / k, 0.0, 1.0);
    return d1 * (1.0 - h) + d2 * h - k * h * (1.0 - h);
}

/**
 * @brief 平滑交集（带混合过渡）
 *
 * 与平滑并集类似，但应用于交集操作。
 *
 * **数学公式:**
 * h = clamp(0.5 - 0.5·(d1-d2)/k, 0, 1)
 * f = d1·(1-h) + d2·h + k·h·(1-h)
 *
 * @param d1 第一个对象的 SDF 值
 * @param d2 第二个对象的 SDF 值
 * @param k 混合强度参数
 * @return 平滑混合后的 SDF 值
 *
 * @see op_intersection 无平滑的标准交集
 * @see op_smooth_union 平滑并集
 */
[[nodiscard]] inline double op_smooth_intersection(double d1, double d2, double k) {
    if (k <= 0.0) return std::max(d1, d2);
    double h = std::clamp(0.5 - 0.5 * (d1 - d2) / k, 0.0, 1.0);
    return d1 * (1.0 - h) + d2 * h + k * h * (1.0 - h);
}

/**
 * @brief 平滑差集（带混合过渡）
 *
 * 与平滑并集类似，但应用于差集操作。
 *
 * **数学公式:**
 * h = clamp(0.5 - 0.5·(d2+d1)/k, 0, 1)
 * f = d1·(1-h) + (-d2)·h + k·h·(1-h)
 *
 * @param d1 被减对象 A 的 SDF 值
 * @param d2 减去对象 B 的 SDF 值
 * @param k 混合强度参数
 * @return 平滑混合后的 SDF 值
 *
 * @see op_difference 无平滑的标准差集
 * @see op_smooth_union 平滑并集
 */
[[nodiscard]] inline double op_smooth_difference(double d1, double d2, double k) {
    if (k <= 0.0) return std::max(d1, -d2);
    double h = std::clamp(0.5 - 0.5 * (d2 + d1) / k, 0.0, 1.0);
    return d1 * (1.0 - h) + (-d2) * h + k * h * (1.0 - h);
}

// ── Modifiers (operate on a distance value) ──

/**
 * @brief 圆角偏移: 将等值面向外扩展 r 单位
 *
 * **数学公式:** f = d - r
 *
 * 相当于对隐式曲面做正向膨胀。配合其他运算可实现倒圆角效果。
 * 例如 box SDF + op_round 即等价于 round_box。
 *
 * @param d SDF 值
 * @param r 偏移距离（正值 = 膨胀，负值 = 收缩）
 * @return 偏移后的 SDF 值
 *
 * @code{.cpp}
 * // 给任意形状加 0.1 的圆角
 * double rounded = op_round(my_sdf(p), 0.1);
 * @endcode
 *
 * @see op_onion 壳层修饰
 * @see round_box 特例化的圆角立方体
 */
[[nodiscard]] inline double op_round(double d, double r) {
    return d - r;
}

/**
 * @brief 洋葱皮/壳层修饰: 取绝对值后减去厚度
 *
 * **数学公式:** f = |d| - thickness
 *
 * 将 SDF 的等值面变为一个厚度为 2·thickness 的壳层。
 * 原始表面变成壳层的内外边界。
 *
 * @param d SDF 值
 * @param thickness 壳层半厚度，必须 > 0
 * @return 壳层的 SDF 值
 *
 * @code{.cpp}
 * // 0.05 厚的球壳
 * double shell = op_onion(sphere(p, 1.0), 0.05);
 * @endcode
 *
 * @see op_round 圆角偏移
 */
[[nodiscard]] inline double op_onion(double d, double thickness) {
    return std::abs(d) - thickness;
}

// ── Domain deformation operators (transform point space) ──

/**
 * @brief 无限重复（周期复制）
 *
 * 将空间映射到以 cell 为周期的单位晶格中。
 * 每个晶格是边长为 cell 的轴对齐长方体。
 *
 * **变换:** p' = p - cell·round(p / cell)
 *
 * 配合任何 SDF 使用，可生成无限周期阵列。
 *
 * @param p 原始采样点
 * @param cell 晶格尺寸 (cx, cy, cz)，各分量必须 > 0 才生效
 * @return 映射到中心晶格内的点
 *
 * @code{.cpp}
 * // 在空间中无限重复的球体阵列（间距 2）
 * Point3D q = op_repeat(p, Point3D(2, 2, 2));
 * double d = sphere(q, 0.5);
 * @endcode
 *
 * @note 边界处 SDF 可能不连续，需要确保各晶格间映射是连续的。
 */
[[nodiscard]] inline Point3D op_repeat(const Point3D& p, const Point3D& cell) {
    auto wrap = [](double v, double c) {
        return c > 0.0 ? v - c * std::round(v / c) : v;
    };
    return Point3D(wrap(p.x(), cell.x()),
                   wrap(p.y(), cell.y()),
                   wrap(p.z(), cell.z()));
}

/**
 * @brief 镜像对称（X=0 平面），可选偏移
 *
 * 将空间关于 x = offset 平面取镜像对称。
 *
 * **变换:** p'_x = offset + |p_x - offset|
 *
 * @param p 原始采样点
 * @param offset 镜像平面沿 X 轴的偏移量（默认 0 = X=0 平面）
 * @return 映射到镜像半空间的点
 *
 * @code{.cpp}
 * // 关于 X=1 平面对称
 * Point3D q = op_mirror_x(p, 1.0);
 * @endcode
 *
 * @see op_mirror_y YZ 平面镜像
 * @see op_mirror_z XY 平面镜像
 */
[[nodiscard]] inline Point3D op_mirror_x(const Point3D& p, double offset = 0.0) {
    return Point3D(offset + std::abs(p.x() - offset), p.y(), p.z());
}

/**
 * @brief 镜像对称（Y=0 平面）
 *
 * 将空间关于 XZ 平面取镜像对称。
 *
 * **变换:** p'_y = |p_y|
 *
 * @param p 原始采样点
 * @return 映射到上半空间的点
 *
 * @code{.cpp}
 * // 关于 Y=0 对称（上下对称）
 * Point3D q = op_mirror_y(p);
 * double d = sphere(q - Point3D(0, 1.0, 0), 0.5);  // 上下各一个球
 * @endcode
 */
[[nodiscard]] inline Point3D op_mirror_y(const Point3D& p) {
    return Point3D(p.x(), std::abs(p.y()), p.z());
}

/**
 * @brief 镜像对称（Z=0 平面）
 *
 * 将空间关于 XY 平面取镜像对称。
 *
 * **变换:** p'_z = |p_z|
 *
 * @param p 原始采样点
 * @return 映射到前半空间的点
 */
[[nodiscard]] inline Point3D op_mirror_z(const Point3D& p) {
    return Point3D(p.x(), p.y(), std::abs(p.z()));
}

/**
 * @brief 空间平移（逆变换）
 *
 * 将采样点向反方向平移，使后续 SDF 求值表现出正向移动效果。
 *
 * **变换:** p' = p - offset
 *
 * 例如 offset = (2, 0, 0) 会使 SDF 在空间中表现为向右移动 2 个单位。
 *
 * @param p 原始采样点
 * @param offset 平移量（对象正向移动的方向）
 * @return 逆平移后的采样点
 *
 * @code{.cpp}
 * // 球体向右平移 2 个单位
 * Point3D q = op_translate(p, Point3D(2, 0, 0));
 * double d = sphere(q, 1.0);
 * @endcode
 *
 * @see op_rotate 旋转变换
 * @see op_scale 缩放变换
 */
[[nodiscard]] inline Point3D op_translate(const Point3D& p, const Point3D& offset) {
    return p - offset;
}

/**
 * @brief 绕 Y 轴旋转（逆变换）
 *
 * 将采样点绕 Y 轴旋转 -angle_rad，使后续 SDF 表现出正向旋转效果。
 *
 * **变换矩阵:**
 * ```
 * p'_x =  p_x·cos(θ) + p_z·sin(θ)
 * p'_y =  p_y
 * p'_z = -p_x·sin(θ) + p_z·cos(θ)
 * ```
 *
 * @param p 原始采样点
 * @param angle_rad 旋转角度（弧度），正值 = 逆时针（从上往下看）
 * @return 逆旋转后的采样点
 *
 * @code{.cpp}
 * // 将盒子绕 Y 轴旋转 45°
 * Point3D q = op_rotate(p, M_PI / 4);
 * double d = box(q, Point3D(1, 1, 1));
 * @endcode
 *
 * @see op_twist 绕 Y 轴的扭曲变换
 * @see op_bend 弯曲变换
 */
[[nodiscard]] inline Point3D op_rotate(const Point3D& p, double angle_rad) {
    double c = std::cos(-angle_rad);
    double s = std::sin(-angle_rad);
    return Point3D(p.x() * c - p.z() * s, p.y(), p.x() * s + p.z() * c);
}

/**
 * @brief 非均匀缩放（逆变换）
 *
 * 将采样点各分量分别除以缩放因子。
 *
 * **变换:** p'_i = p_i / s_i
 *
 * @param p 原始采样点
 * @param s 各轴缩放因子 (sx, sy, sz)，均必须 ≠ 0
 * @return 缩放后的采样点
 *
 * @warning 非均匀缩放会破坏 SDF 的距离性质（梯度不再是单位长度）。
 *          只有均匀缩放（sx=sy=sz）才能保持距离场精确性。
 *
 * @code{.cpp}
 * // 均匀缩小为一半
 * Point3D q = op_scale(p, Point3D(2, 2, 2));
 * double d = sphere(q, 1.0);  // 实际半径变为 0.5
 * @endcode
 *
 * @see op_translate 平移（等距变换）
 * @see op_rotate 旋转（等距变换）
 */
[[nodiscard]] inline Point3D op_scale(const Point3D& p, const Point3D& s) {
    return Point3D(p.x() / s.x(), p.y() / s.y(), p.z() / s.z());
}

/**
 * @brief 绕 Y 轴的扭曲变换
 *
 * 旋转角度随 Y 坐标线性变化: angle = amount · p_y。
 *
 * **变换:**
 * ```
 * p'_x =  p_x·cos(amount·p_y) - p_z·sin(amount·p_y)
 * p'_y =  p_y
 * p'_z =  p_x·sin(amount·p_y) + p_z·cos(amount·p_y)
 * ```
 *
 * @param p 原始采样点
 * @param amount 单位高度上的旋转量（弧度/单位长度）
 * @return 扭曲后的采样点
 *
 * @code{.cpp}
 * // 在上方 2 个单位处旋转 90 度的扭曲
 * Point3D q = op_twist(p, M_PI / 4);
 * double d = box(q, Point3D(0.5, 2, 0.5));
 * @endcode
 *
 * @see op_rotate 恒定角度旋转
 * @see op_bend 弯曲变换
 */
[[nodiscard]] inline Point3D op_twist(const Point3D& p, double amount) {
    double c = std::cos(amount * p.y());
    double s = std::sin(amount * p.y());
    return Point3D(p.x() * c - p.z() * s, p.y(), p.x() * s + p.z() * c);
}

/**
 * @brief 弯曲变换（绕 Z 轴弯曲 XZ 平面）
 *
 * 将空间沿 X 轴弯曲，弯曲量由参数 k 控制。
 *
 * **变换:**
 * ```
 * p'_x =  cos(k·p_x)·p_x - sin(k·p_x)·p_y
 * p'_y =  sin(k·p_x)·p_x + cos(k·p_x)·p_y
 * p'_z =  p_z
 * ```
 *
 * @param p 原始采样点
 * @param k 弯曲曲率参数（k 越大弯曲越剧烈）
 * @return 弯曲后的采样点
 *
 * @code{.cpp}
 * // 轻度弯曲的长方体
 * Point3D q = op_bend(p, 0.2);
 * double d = box(q, Point3D(3, 0.3, 0.3));
 * @endcode
 *
 * @see op_cheap_bend 简化版弯曲（仅绕 X 轴）
 * @see op_twist 扭曲变换
 */
[[nodiscard]] inline Point3D op_bend(const Point3D& p, double k) {
    double c = std::cos(k * p.x());
    double s = std::sin(k * p.x());
    return Point3D(c * p.x() - s * p.y(), s * p.x() + c * p.y(), p.z());
}

/**
 * @brief 拉伸/压扁域变形（消除拉伸方向上的内部体积）
 *
 * 将点沿各主轴方向向原点压缩，消除 |p_i| ≤ h_i 区域的内部信息。
 * 这在组合建模中用于从对象中挖去沿坐标轴的孔道。
 *
 * **变换:**
 * 对每个分量 i：
 * - 若 |p_i| > h_i: p'_i = p_i - sign(p_i)·h_i
 * - 否则: p'_i = 0
 *
 * @param p 原始采样点
 * @param h 各轴的拉伸阈值 (hx, hy, hz)，均 ≥ 0
 * @return 变形后的采样点
 *
 * @code{.cpp}
 * // 将球体沿 X 轴拉伸，挖去|x|≤1 的区域
 * Point3D q = op_elongate(p, Point3D(1, 0, 0));
 * double d = sphere(q, 1.5);
 * @endcode
 *
 * @see op_repeat 周期重复
 */
[[nodiscard]] inline Point3D op_elongate(const Point3D& p, const Point3D& h) {
    Point3D q = p.cwiseAbs() - h;
    return Point3D(
        q.x() > 0.0 ? p.x() - std::copysign(h.x(), p.x()) : 0.0,
        q.y() > 0.0 ? p.y() - std::copysign(h.y(), p.y()) : 0.0,
        q.z() > 0.0 ? p.z() - std::copysign(h.z(), p.z()) : 0.0
    );
}

/**
 * @brief 简化弯曲变换（绕 X 轴弯曲 YZ 平面）
 *
 * 比 op_bend 更快但精度较低的弯曲，仅在 YZ 方向产生偏移。
 *
 * **变换:**
 * ```
 * p'_x = p_x
 * p'_y = p_y·cos(k·p_x) - p_z·sin(k·p_x)
 * p'_z = p_y·sin(k·p_x) + p_z·cos(k·p_x)
 * ```
 *
 * @param p 原始采样点
 * @param k 弯曲曲率参数
 * @return 弯曲后的采样点
 *
 * @note 使用场景：不需要精确物理弯曲，仅需视觉效果时优先用此函数以节省计算。
 *
 * @see op_bend 完整弯曲变换
 */
[[nodiscard]] inline Point3D op_cheap_bend(const Point3D& p, double k) {
    double c = std::cos(k * p.x());
    double s = std::sin(k * p.x());
    return Point3D(p.x(), p.y() * c - p.z() * s, p.y() * s + p.z() * c);
}

/**
 * @brief 正弦波位移扰动
 *
 * 在距离场上叠加一个三维正弦波的微小扰动。
 *
 * **数学公式:**
 * f = d + amplitude · ∏_i sin(p_i · frequency)
 *
 * 三个方向的 sin 函数乘积产生体纹理效果。
 *
 * @param d 原始 SDF 值
 * @param p 采样点（用于扰动计算）
 * @param amplitude 扰动幅值
 * @param frequency 空间频率
 * @return 扰动后的 SDF 值
 *
 * @warning 扰动能轻易破坏距离场的 Lipschitz 性质（|∇f| ≤ 1），
 *          导致 ray-marching 效率下降。需控制 amplitude/frequency 在合理范围内。
 *
 * @code{.cpp}
 * // 给球体添加微小波纹
 * double d = op_displace(sphere(p, 1.0), p, 0.05, 4.0);
 * @endcode
 */
[[nodiscard]] inline double op_displace(double d, const Point3D& p,
                                         double amplitude, double frequency) {
    double noise = std::sin(p.x() * frequency) *
                   std::sin(p.y() * frequency) *
                   std::sin(p.z() * frequency);
    return d + amplitude * noise;
}

} // namespace vde::sdf