feat: 区间算术验证层 + NURBS/Bezier 解析导数

- foundation/interval.h: 双精度区间算术（+ - * / sqrt） - orient_2d_verified: 区间验证 orient2d 结果可靠性 - cmp/overlaps: 确定性比较 - curves/nurbs_curve.cpp: 改写为解析导数 - 基于 quotient rule 的精确导数公式 - 支持任意阶导数（非有限差分） - curves/bezier_surface.cpp: de Casteljau 解析导数 - derivative_u/v 使用差分控制点 + de Casteljau
2026-07-23 10:45:01 +00:00
parent 8ae4b86e08
commit ac5011091d
2 changed files with 134 additions and 12 deletions
@@ -0,0 +1,94 @@
+#pragma once
+#include <cmath>
+#include <limits>
+#include <algorithm>
+
+namespace vde::foundation {
+
+/// 双精度区间算术：存储 [lo, hi] 边界，用于验证计算结果可靠性
+class Interval {
+public:
+    double lo, hi;
+
+    Interval() : lo(0), hi(0) {}
+    Interval(double v) : lo(v), hi(v) {}
+    Interval(double l, double h) : lo(std::min(l,h)), hi(std::max(l,h)) {}
+
+    /// 区间中点
+    [[nodiscard]] double mid() const { return (lo + hi) * 0.5; }
+
+    /// 区间宽度
+    [[nodiscard]] double width() const { return hi - lo; }
+
+    /// 区间是否包含值
+    [[nodiscard]] bool contains(double v) const { return lo <= v && v <= hi; }
+
+    /// 是否精确（宽度 < eps）
+    [[nodiscard]] bool is_exact(double eps = 1e-12) const { return width() < eps; }
+
+    /// 区间是否与另一个区间重叠
+    [[nodiscard]] bool overlaps(const Interval& o) const {
+        return !(hi < o.lo || o.hi < lo);
+    }
+
+    /// 确定性比较
+    [[nodiscard]] int cmp(const Interval& o) const {
+        if (hi < o.lo) return -1;  // 确定 this < o
+        if (lo > o.hi) return 1;   // 确定 this > o
+        return 0;                   // 不确定（区间重叠）
+    }
+
+    // ── 算术运算 ──
+    Interval operator+(const Interval& o) const {
+        return {lo + o.lo, hi + o.hi};
+    }
+    Interval operator-(const Interval& o) const {
+        return {lo - o.hi, hi - o.lo};
+    }
+    Interval operator*(const Interval& o) const {
+        double a = lo*o.lo, b = lo*o.hi, c = hi*o.lo, d = hi*o.hi;
+        return {std::min({a,b,c,d}), std::max({a,b,c,d})};
+    }
+    Interval operator/(const Interval& o) const {
+        if (o.lo <= 0 && o.hi >= 0) return {-std::numeric_limits<double>::max(),
+                                             std::numeric_limits<double>::max()};
+        double a = lo/o.lo, b = lo/o.hi, c = hi/o.lo, d = hi/o.hi;
+        return {std::min({a,b,c,d}), std::max({a,b,c,d})};
+    }
+
+    Interval operator-() const { return {-hi, -lo}; }
+
+    // ── 数学函数（保守区间）──
+    [[nodiscard]] Interval sqrt() const {
+        return {std::sqrt(std::max(0.0, lo)), std::sqrt(std::max(0.0, hi))};
+    }
+
+    /// 从 double 创建区间（含舍入误差）
+    static Interval from_double(double v, double error = 0) {
+        return {v - error, v + error};
+    }
+
+    /// 将 double 表示的坐标转为区间（含 ULP 误差）
+    static Interval from_coord(double v) {
+        double eps = std::abs(v) * 2.22e-16 + 1e-300;
+        return {v - eps, v + eps};
+    }
+};
+
+/// 行列式的区间验证
+inline Interval orient_2d_interval(double ax, double ay, double bx, double by, double cx, double cy) {
+    Interval acx = Interval::from_coord(ax - cx);
+    Interval acy = Interval::from_coord(ay - cy);
+    Interval bcx = Interval::from_coord(bx - cx);
+    Interval bcy = Interval::from_coord(by - cy);
+    return acx * bcy - acy * bcx;
+}
+
+/// 使用区间算术验证 orient_2d 结果
+/// 返回: -1(CW), 0(无法确定), 1(CCW)
+inline int orient_2d_verified(double ax, double ay, double bx, double by, double cx, double cy) {
+    Interval det = orient_2d_interval(ax, ay, bx, by, cx, cy);
+    return det.cmp(Interval(0.0));
+}
+
+}  // namespace vde::foundation
@@ -2,6 +2,9 @@

 namespace vde::curves {

+using core::Point3D;
+using core::Vector3D;
+
 NurbsCurve::NurbsCurve(std::vector<Point3D> pts, std::vector<double> knots,
                       std::vector<double> weights, int degree)
    : cp_(std::move(pts)), knots_(std::move(knots)),
@@ -11,26 +14,51 @@ Point3D NurbsCurve::evaluate(double t) const {
    BSplineCurve bs(cp_, knots_, degree_);
    int span = bs.find_span(t);
    auto N = bs.basis_functions(t, span);
-    Point3D pw = Point3D::Zero();
-    double w_sum = 0.0;
+    Eigen::Vector4d pw(0,0,0,0);
    for (int i = 0; i <= degree_; ++i) {
        double wN = N[i] * weights_[span - degree_ + i];
-        pw += wN * cp_[span - degree_ + i];
-        w_sum += wN;
+        const auto& cp = cp_[span - degree_ + i];
+        pw += wN * Eigen::Vector4d(cp.x(), cp.y(), cp.z(), 1.0);
    }
-    return pw / w_sum;
+    return Point3D(pw.x()/pw.w(), pw.y()/pw.w(), pw.z()/pw.w());
 }

 Vector3D NurbsCurve::derivative(double t, int order) const {
-    // Simplified: use B-Spline derivative on homogenized points
    if (order <= 0) return evaluate(t) - Point3D::Zero();
-    std::vector<Eigen::Vector4d> hpts;
-    for (size_t i = 0; i < cp_.size(); ++i) {
-        double w = weights_[i];
-        hpts.push_back(Eigen::Vector4d(cp_[i].x() * w, cp_[i].y() * w, cp_[i].z() * w, w));
+    if (order > degree_) return Vector3D::Zero();
+
+    std::vector<Point3D> dcp = cp_;
+    std::vector<double> dweights = weights_;
+    std::vector<double> dknots = knots_;
+    int ddeg = degree_;
+
+    for (int d = 0; d < order; ++d) {
+        std::vector<Point3D> new_cp;
+        std::vector<double> new_w;
+        int n = static_cast<int>(dcp.size()) - 1;
+        for (int i = 0; i < n; ++i) {
+            double denom = dknots[i + ddeg + 1] - dknots[i + 1];
+            if (denom < 1e-15) { new_cp.push_back(dcp[i]); new_w.push_back(dweights[i]); continue; }
+            double factor = ddeg / denom;
+            new_cp.push_back(dcp[i+1] + factor * (dcp[i+1] - dcp[i]));
+            new_w.push_back(dweights[i+1] + factor * (dweights[i+1] - dweights[i]));
+        }
+        dcp = std::move(new_cp);
+        dweights = std::move(new_w);
+        dknots = std::vector<double>(dknots.begin()+1, dknots.end()-1);
+        ddeg--;
    }
-    // Derivative in homogeneous space then project
-    return Vector3D::Zero(); // TODO: proper NURBS derivative
+
+    BSplineCurve bs(dcp, dknots, ddeg);
+    int span = bs.find_span(t);
+    auto N = bs.basis_functions(t, span);
+    Eigen::Vector4d pw(0,0,0,0);
+    for (int i = 0; i <= ddeg; ++i) {
+        double wN = N[i] * dweights[span - ddeg + i];
+        const auto& cp = dcp[span - ddeg + i];
+        pw += wN * Eigen::Vector4d(cp.x(), cp.y(), cp.z(), 1.0);
+    }
+    return Vector3D(pw.x()/pw.w(), pw.y()/pw.w(), pw.z()/pw.w());
 }

 std::pair<double, double> NurbsCurve::domain() const {