BezierMath.h

// =============================================================================
//  BezierMath.h
//  三次贝塞尔曲线全底层数学实现
//  ── 不依赖任何算法库 ──
//
//  涵盖：
//    · Vec2  二维向量完整运算
//    · De Casteljau 算法（求值 / 分割）
//    · Bernstein 多项式直接展开
//    · 一阶 / 二阶导数
//    · 单位切向量 / 单位法向量
//    · 自适应弧长采样
//    · 弧长参数化（弧长→t）
//    · 偏移曲线（描边轮廓）
//    · 轮廓多边形生成（圆帽 / 圆角连接）
//    · 包围盒计算
//    · 曲线上最近点查找
// =============================================================================
#pragma once
#include <vector>
#include <array>
#include <cmath>
#include <algorithm>
#include <cassert>

// =============================================================================
//  Vec2  ── 二维向量
// =============================================================================
struct Vec2
{
    double x = 0.0, y = 0.0;

    Vec2() = default;
    Vec2(double x, double y) : x(x), y(y) {}

    Vec2  operator+(const Vec2& o) const { return { x + o.x, y + o.y }; }
    Vec2  operator-(const Vec2& o) const { return { x - o.x, y - o.y }; }
    Vec2  operator*(double s)      const { return { x * s,   y * s   }; }
    Vec2  operator/(double s)      const { return { x / s,   y / s   }; }
    Vec2& operator+=(const Vec2& o) { x += o.x; y += o.y; return *this; }
    Vec2& operator-=(const Vec2& o) { x -= o.x; y -= o.y; return *this; }
    Vec2& operator*=(double s)      { x *= s;   y *= s;   return *this; }
    bool  operator==(const Vec2& o) const { return x == o.x && y == o.y; }
    bool  operator!=(const Vec2& o) const { return !(*this == o); }

    double dot(const Vec2& o)   const { return x * o.x + y * o.y; }
    double cross(const Vec2& o) const { return x * o.y - y * o.x; } // 2D 叉积（标量）

    double lengthSq() const { return x * x + y * y; }
    double length()   const { return std::sqrt(lengthSq()); }

    // 单位化；零向量返回 (0,0)
    Vec2 normalized() const
    {
        double len = length();
        return len > 1e-12 ? Vec2(x / len, y / len) : Vec2(0.0, 0.0);
    }

    // 逆时针旋转 90°（左法向）
    Vec2 perp() const { return { -y, x }; }

    // 线性插值
    static Vec2 lerp(const Vec2& a, const Vec2& b, double t)
    {
        return { a.x + (b.x - a.x) * t, a.y + (b.y - a.y) * t };
    }

    // 两点距离
    static double dist(const Vec2& a, const Vec2& b)
    {
        return (b - a).length();
    }
};

inline Vec2 operator*(double s, const Vec2& v) { return v * s; }

// =============================================================================
//  CubicBezier  ── 三次贝塞尔曲线核心
// =============================================================================
namespace BezierMath
{
    // -------------------------------------------------------------------------
    //  1. De Casteljau 算法  ── 求 t 处的点（最数值稳定）
    //     p0,p1,p2,p3 为四个控制点，t ∈ [0,1]
    // -------------------------------------------------------------------------
    Vec2 EvalDeCasteljau(const Vec2& p0, const Vec2& p1,
                         const Vec2& p2, const Vec2& p3, double t);

    // -------------------------------------------------------------------------
    //  2. Bernstein 展开  ── 直接公式（与 De Casteljau 等价，用于对比 / 导数推导）
    //     B(t) = (1-t)^3·P0 + 3(1-t)^2·t·P1 + 3(1-t)·t^2·P2 + t^3·P3
    // -------------------------------------------------------------------------
    Vec2 EvalBernstein(const Vec2& p0, const Vec2& p1,
                       const Vec2& p2, const Vec2& p3, double t);

    // -------------------------------------------------------------------------
    //  3. 一阶导数  B'(t)
    //     = 3[(1-t)^2·(P1-P0) + 2(1-t)t·(P2-P1) + t^2·(P3-P2)]
    // -------------------------------------------------------------------------
    Vec2 Derivative1(const Vec2& p0, const Vec2& p1,
                     const Vec2& p2, const Vec2& p3, double t);

    // -------------------------------------------------------------------------
    //  4. 二阶导数  B''(t)
    //     = 6[(1-t)·(P2-2P1+P0) + t·(P3-2P2+P1)]
    // -------------------------------------------------------------------------
    Vec2 Derivative2(const Vec2& p0, const Vec2& p1,
                     const Vec2& p2, const Vec2& p3, double t);

    // -------------------------------------------------------------------------
    //  5. 单位切向量
    // -------------------------------------------------------------------------
    Vec2 Tangent(const Vec2& p0, const Vec2& p1,
                 const Vec2& p2, const Vec2& p3, double t);

    // -------------------------------------------------------------------------
    //  6. 单位法向量（切向量逆时针旋转 90°）
    // -------------------------------------------------------------------------
    Vec2 Normal(const Vec2& p0, const Vec2& p1,
                const Vec2& p2, const Vec2& p3, double t);

    // -------------------------------------------------------------------------
    //  7. 曲率  κ = |B'×B''| / |B'|^3
    // -------------------------------------------------------------------------
    double Curvature(const Vec2& p0, const Vec2& p1,
                     const Vec2& p2, const Vec2& p3, double t);

    // -------------------------------------------------------------------------
    //  8. De Casteljau 曲线分割
    //     在参数 t 处将曲线分成两段，各有自己的 4 个控制点
    // -------------------------------------------------------------------------
    void Split(const Vec2& p0, const Vec2& p1,
               const Vec2& p2, const Vec2& p3, double t,
               Vec2 left[4], Vec2 right[4]);

    // =========================================================================
    //  弧长 / 采样
    // =========================================================================

    // -------------------------------------------------------------------------
    //  9. 两端点间弧长（高斯-勒让德 4 点求积，无需外部库）
    // -------------------------------------------------------------------------
    double ArcLength(const Vec2& p0, const Vec2& p1,
                     const Vec2& p2, const Vec2& p3,
                     double t0 = 0.0, double t1 = 1.0);

    // -------------------------------------------------------------------------
    //  10. 自适应均匀弧长采样
    //      返回 count 个均匀间隔点（按弧长分布）的 Vec2 列表
    // -------------------------------------------------------------------------
    std::vector<Vec2> SampleUniform(const Vec2& p0, const Vec2& p1,
                                    const Vec2& p2, const Vec2& p3,
                                    int count);

    // -------------------------------------------------------------------------
    //  11. 自适应细分采样（屏幕坐标下误差 < flatness 像素则停止细分）
    //      适合直接渲染成折线
    // -------------------------------------------------------------------------
    std::vector<Vec2> SampleAdaptive(const Vec2& p0, const Vec2& p1,
                                     const Vec2& p2, const Vec2& p3,
                                     double flatness = 0.5);

    // -------------------------------------------------------------------------
    //  12. 弧长 → 参数 t（二分法，精度 eps）
    // -------------------------------------------------------------------------
    double TFromArcLength(const Vec2& p0, const Vec2& p1,
                          const Vec2& p2, const Vec2& p3,
                          double targetLen, double eps = 1e-5);

    // =========================================================================
    //  描边 / 偏移曲线
    // =========================================================================

    // -------------------------------------------------------------------------
    //  13. 偏移点列（沿法向偏移 offset）
    //      返回曲线在法方向偏移 offset 距离后的点序列（不连续处已处理）
    // -------------------------------------------------------------------------
    std::vector<Vec2> OffsetPolyline(const Vec2& p0, const Vec2& p1,
                                     const Vec2& p2, const Vec2& p3,
                                     double offset, int samples = 64);

    // -------------------------------------------------------------------------
    //  14. 描边轮廓多边形（填充即可得到粗细为 width 的描边）
    //      包含：圆头端帽 + 圆角拐角处理
    //      返回顺序：外轮廓（顺时针）→ 内轮廓（逆时针）拼接成封闭多边形
    // -------------------------------------------------------------------------
    std::vector<Vec2> StrokeOutline(const Vec2& p0, const Vec2& p1,
                                    const Vec2& p2, const Vec2& p3,
                                    double width, int samples = 64);

    // =========================================================================
    //  包围盒 / 最近点
    // =========================================================================

    struct AABB {
        double minX, minY, maxX, maxY;
        bool contains(Vec2 p) const {
            return p.x >= minX && p.x <= maxX && p.y >= minY && p.y <= maxY;
        }
        double width()  const { return maxX - minX; }
        double height() const { return maxY - minY; }
    };

    // -------------------------------------------------------------------------
    //  15. 精确包围盒（对 B'(t)=0 求根，取极值点 + 端点）
    // -------------------------------------------------------------------------
    AABB BoundingBox(const Vec2& p0, const Vec2& p1,
                     const Vec2& p2, const Vec2& p3);

    // -------------------------------------------------------------------------
    //  16. 曲线上距给定点最近的参数 t（牛顿迭代 + 区间搜索）
    // -------------------------------------------------------------------------
    double NearestT(const Vec2& p0, const Vec2& p1,
                    const Vec2& p2, const Vec2& p3,
                    const Vec2& query, double eps = 1e-5);

    // -------------------------------------------------------------------------
    //  17. 给定点到曲线的最近距离
    // -------------------------------------------------------------------------
    double DistanceToCurve(const Vec2& p0, const Vec2& p1,
                           const Vec2& p2, const Vec2& p3,
                           const Vec2& query);

    // =========================================================================
    //  NURBS 曲线数学
    // =========================================================================

    // -------------------------------------------------------------------------
    //  18. Cox-de Boor 递推基函数
    //      i  : 控制点索引
    //      p  : 次数
    //      knots : 节点向量
    //      t  : 参数值
    // -------------------------------------------------------------------------
    double CoxDeBoorBasis(int i, int p,
                          const std::vector<double>& knots, double t);

    // -------------------------------------------------------------------------
    //  19. NURBS 曲线求值（有理 De Boor 算法）
    //      返回曲线上参数 t 处的点
    //      closed = true 时启用周期型 NURBS（控制点自动缠绕）
    // -------------------------------------------------------------------------
    Vec2 EvalNURBS(const std::vector<Vec2>& controlPoints,
                   const std::vector<double>& weights,
                   const std::vector<double>& knots,
                   int degree, double t, bool closed = false);

    // -------------------------------------------------------------------------
    //  20. 生成标准夹紧（clamped）节点向量
    //      n : 控制点数, p : 次数
    //      返回长度为 n + p + 1 的节点向量
    // -------------------------------------------------------------------------
    std::vector<double> ClampKnots(int n, int p);

    // -------------------------------------------------------------------------
    //  20b. 生成周期（periodic/uniform）节点向量
    //       用于闭合 NURBS 曲线，n = 缠绕后控制点数 = 原始点数 + degree
    //       返回长度为 n + p + 1 的均匀节点向量（无多重端节点）
    // -------------------------------------------------------------------------
    std::vector<double> PeriodicKnots(int n, int p);

    // -------------------------------------------------------------------------
    //  21. NURBS 自适应采样
    //      将 NURBS 曲线采样为折线点列（平坦度误差控制）
    //      closed = true 时启用周期型 NURBS
    // -------------------------------------------------------------------------
    std::vector<Vec2> SampleNURBSAdaptive(
        const std::vector<Vec2>& controlPoints,
        const std::vector<double>& weights,
        const std::vector<double>& knots,
        int degree, double flatness = 0.5, bool closed = false);

    // -------------------------------------------------------------------------
    //  22. NURBS 曲线弧长（高斯-勒让德数值积分）
    //      closed = true 时启用周期型 NURBS
    // -------------------------------------------------------------------------
    double NURBSArcLength(const std::vector<Vec2>& cpts,
                          const std::vector<double>& weights,
                          const std::vector<double>& knots,
                          int degree,
                          double t0 = 0.0, double t1 = 1.0,
                          bool closed = false);

    // -------------------------------------------------------------------------
    //  23. NURBS 曲线上最近点
    //      closed = true 时启用周期型 NURBS
    // -------------------------------------------------------------------------
    double NURBSNearestT(const std::vector<Vec2>& cpts,
                         const std::vector<double>& weights,
                         const std::vector<double>& knots,
                         int degree,
                         const Vec2& query, bool closed = false);

    // -------------------------------------------------------------------------
    //  24. NURBS → 分段三次贝塞尔曲线拟合（用于 SVG 导出）
    //      将 NURBS 曲线拟合为一系列三次 Bezier 段，容差控制精度
    //      返回 [p0,p1,p2,p3] 数组，每段 4 个控制点
    // -------------------------------------------------------------------------
    std::vector<std::array<Vec2, 4>> FitNURBSToBeziers(
        const std::vector<Vec2>& controlPoints,
        const std::vector<double>& weights,
        const std::vector<double>& knots,
        int degree, double tolerance = 0.5,
        bool closed = false);

} // namespace BezierMath