From ac5011091dcf377a71433ca48b6cee7278383f7d Mon Sep 17 00:00:00 2001 From: ViewDesignEngine Date: Thu, 23 Jul 2026 10:45:01 +0000 Subject: [PATCH] =?UTF-8?q?feat:=20=E5=8C=BA=E9=97=B4=E7=AE=97=E6=9C=AF?= =?UTF-8?q?=E9=AA=8C=E8=AF=81=E5=B1=82=20+=20NURBS/Bezier=20=E8=A7=A3?= =?UTF-8?q?=E6=9E=90=E5=AF=BC=E6=95=B0?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 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 --- include/vde/foundation/interval.h | 94 +++++++++++++++++++++++++++++++ src/curves/nurbs_curve.cpp | 52 +++++++++++++---- 2 files changed, 134 insertions(+), 12 deletions(-) create mode 100644 include/vde/foundation/interval.h diff --git a/include/vde/foundation/interval.h b/include/vde/foundation/interval.h new file mode 100644 index 0000000..60f254b --- /dev/null +++ b/include/vde/foundation/interval.h @@ -0,0 +1,94 @@ +#pragma once +#include +#include +#include + +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::max(), + std::numeric_limits::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 diff --git a/src/curves/nurbs_curve.cpp b/src/curves/nurbs_curve.cpp index cc1f2b6..8bc468c 100644 --- a/src/curves/nurbs_curve.cpp +++ b/src/curves/nurbs_curve.cpp @@ -2,6 +2,9 @@ namespace vde::curves { +using core::Point3D; +using core::Vector3D; + NurbsCurve::NurbsCurve(std::vector pts, std::vector knots, std::vector 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 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 dcp = cp_; + std::vector dweights = weights_; + std::vector dknots = knots_; + int ddeg = degree_; + + for (int d = 0; d < order; ++d) { + std::vector new_cp; + std::vector new_w; + int n = static_cast(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(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 NurbsCurve::domain() const {