diff --git a/include/vde/core/aabb.h b/include/vde/core/aabb.h index f838f21..c49737f 100644 --- a/include/vde/core/aabb.h +++ b/include/vde/core/aabb.h @@ -5,50 +5,131 @@ namespace vde::core { +/** + * @brief 轴对齐包围盒（AABB） + * + * 用最小/最大两个角点表示的三维轴对齐包围盒。 + * 支持增量扩展、包含/相交测试、几何属性查询。 + * + * @tparam T 标量类型（通常为 double 或 float） + * + * @note 初始状态为空（min = +∞, max = -∞），需要通过 expand() 初始化 + * + * @ingroup core + */ template class AABB { public: + /** + * @brief 构造空包围盒 + * + * min 初始化为 T 最大值，max 初始化为 T 最小值， + * 表示"无内容"的空状态。 + */ AABB() : min_(Point3D::Constant(std::numeric_limits::max())), max_(Point3D::Constant(std::numeric_limits::lowest())) {} + /** + * @brief 从两个角点构造包围盒 + * + * @param min 最小角点 + * @param max 最大角点 + * + * @pre min 的各分量 ≤ max 的各分量 + */ AABB(const Point3D& min, const Point3D& max) : min_(min), max_(max) {} + /** + * @brief 扩展包围盒以包含给定点 + * + * 逐分量取 min/max，保证包围盒包含该点。 + * 若包围盒为空，此调用等效于初始化为包含该点的退化盒。 + * + * @param p 要包含的点 + * + * @code{.cpp} + * AABB3D box; + * for (auto& p : points) box.expand(p); + * @endcode + */ void expand(const Point3D& p) { min_ = min_.cwiseMin(p); max_ = max_.cwiseMax(p); } + /** + * @brief 扩展包围盒以包含另一个 AABB + * + * @param other 另一个包围盒 + */ void expand(const AABB& other) { min_ = min_.cwiseMin(other.min_); max_ = max_.cwiseMax(other.max_); } + /** + * @brief 包围盒中心点 + * @return (min + max) / 2 + */ [[nodiscard]] Point3D center() const { return (min_ + max_) * 0.5; } + + /** + * @brief 包围盒对角线向量（从 min 到 max） + * @return max - min + */ [[nodiscard]] Vector3D extent() const { return max_ - min_; } + + /** + * @brief 包围盒表面积 + * @return 2 * (dx*dy + dy*dz + dz*dx) + */ [[nodiscard]] T surface_area() const { Vector3D e = extent(); return 2.0 * (e.x() * e.y() + e.y() * e.z() + e.z() * e.x()); } + + /** + * @brief 包围盒体积 + * @return dx * dy * dz + */ [[nodiscard]] T volume() const { Vector3D e = extent(); return e.x() * e.y() * e.z(); } + + /** + * @brief 测试点是否在包围盒内部（含边界） + * + * @param p 测试点 + * @return true 若 min ≤ p ≤ max（各分量） + */ [[nodiscard]] bool contains(const Point3D& p) const { return (p.array() >= min_.array()).all() && (p.array() <= max_.array()).all(); } + + /** + * @brief 测试两个包围盒是否相交（含边界接触） + * + * @param other 另一个包围盒 + * @return true 若存在交集 + */ [[nodiscard]] bool intersects(const AABB& other) const { return (min_.array() <= other.max_.array()).all() && (max_.array() >= other.min_.array()).all(); } + /// 获取最小角点 [[nodiscard]] const Point3D& min() const { return min_; } + /// 获取最大角点 [[nodiscard]] const Point3D& max() const { return max_; } private: Point3D min_, max_; }; +/// 双精度三维包围盒 using AABB3D = AABB; +/// 单精度三维包围盒 using AABB3f = AABB; } // namespace vde::core diff --git a/include/vde/core/convex_hull.h b/include/vde/core/convex_hull.h index fd32e5d..c90095a 100644 --- a/include/vde/core/convex_hull.h +++ b/include/vde/core/convex_hull.h @@ -4,11 +4,53 @@ namespace vde::core { -/// 2D convex hull — Graham scan, O(n log n) +/** + * @brief 二维凸包计算（Graham Scan 算法） + * + * 使用 Graham 扫描算法计算点集的凸包，时间复杂度 O(n log n)。 + * 以极角排序点为预处理，然后用栈维护凸包边界。 + * + * @param points 输入点集（至少 3 个不共线点） + * @return 凸包顶点序列（逆时针排列），不含共线的中间点 + * + * @note 自动过滤共线点，返回严格的凸包顶点 + * @note 若所有点共线，返回两个端点 + * + * @code{.cpp} + * std::vector pts = {{0,0}, {1,0}, {1,1}, {0,1}, {0.5,0.5}}; + * auto hull = convex_hull_2d(pts); + * // hull = [(0,0), (1,0), (1,1), (0,1)] + * @endcode + * + * @see convex_hull_3d, Polygon2D::convex_hull_2d + * @ingroup core + */ std::vector convex_hull_2d(const std::vector& points); -/// 3D convex hull — QuickHull -/// Returns triangles of the hull (counter-clockwise when viewed from outside) +/** + * @brief 三维凸包计算（QuickHull 算法） + * + * 使用 QuickHull 算法计算三维点集的凸包。 + * 先构造初始四面体，然后递归地将外部点添加到凸包上。 + * 平均复杂度 O(n log n)，最坏情况 O(n²)。 + * + * @param points 输入点集（至少 4 个不共面点） + * @return 凸包三角形列表，每个三角形三顶点按逆时针排列（从外部观察） + * + * @note 返回的三角形法线方向朝外 + * @note 共面/退化点被自动丢弃 + * + * @code{.cpp} + * std::vector pts = /* 3D 点云 */; + * auto hull_tris = convex_hull_3d(pts); + * for (auto& tri : hull_tris) { + * // tri = {p0, p1, p2} — CCW when viewed from outside + * } + * @endcode + * + * @see convex_hull_2d + * @ingroup core + */ std::vector> convex_hull_3d(const std::vector& points); } // namespace vde::core diff --git a/include/vde/core/distance.h b/include/vde/core/distance.h index deb57fc..3e83dd0 100644 --- a/include/vde/core/distance.h +++ b/include/vde/core/distance.h @@ -8,23 +8,194 @@ namespace vde::core { // ── Point distances ────────────────────────── + +/** + * @brief 两点之间的欧氏距离 + * + * @param a 第一个点 + * @param b 第二个点 + * @return ||b - a|| + * + * @code{.cpp} + * double d = distance(Point3D(0,0,0), Point3D(3,4,0)); // d == 5.0 + * @endcode + * + * @ingroup core + */ double distance(const Point3D& a, const Point3D& b); + +/** + * @brief 点到直线的距离 + * + * 计算点到无限直线的垂直距离。 + * + * @param p 测试点 + * @param line 直线 + * @return 垂直距离 + * + * @ingroup core + */ double distance(const Point3D& p, const Line3D& line); + +/** + * @brief 点到线段的距离 + * + * 若投影点落在线段内部，返回垂直距离；否则返回到最近端点的距离。 + * + * @param p 测试点 + * @param seg 线段 + * @return 最短距离 + * + * @ingroup core + */ double distance(const Point3D& p, const Segment3D& seg); + +/** + * @brief 点到平面的距离 + * + * 返回绝对值（总是非负）。 + * + * @param p 测试点 + * @param plane 平面 + * @return 绝对垂直距离 + * + * @see Plane::signed_distance + * @ingroup core + */ double distance(const Point3D& p, const Plane& plane); + +/** + * @brief 点到三角形的距离 + * + * 若投影点在三角形内部，返回垂直距离；否则返回到最近的边或顶点的距离。 + * + * @param p 测试点 + * @param tri 三角形 + * @return 最短距离 + * + * @ingroup core + */ double distance(const Point3D& p, const Triangle& tri); + +/** + * @brief 点到包围盒的距离 + * + * 若点在包围盒内，返回 0；否则返回到最近面的距离。 + * + * @param p 测试点 + * @param box 包围盒 + * @return 最短距离（内部点返回 0） + * + * @ingroup core + */ double distance(const Point3D& p, const AABB& box); // ── Feature-pair distances ─────────────────── + +/** + * @brief 两条线段之间的最短距离 + * + * 处理平行、相交、最近端点到端点等多种情况。 + * + * @param a 第一条线段 + * @param b 第二条线段 + * @return 最短距离（若相交则为 0） + * + * @ingroup core + */ double distance(const Segment3D& a, const Segment3D& b); + +/** + * @brief 两条直线之间的最短距离 + * + * 若相交则为 0；若平行则返回任意垂线段长度；若异面则返回公垂线长度。 + * + * @param a 第一条直线 + * @param b 第二条直线 + * @return 最短距离 + * + * @ingroup core + */ double distance(const Line3D& a, const Line3D& b); + +/** + * @brief 两个包围盒之间的最短距离 + * + * 若相交则返回 0。 + * + * @param a 第一个包围盒 + * @param b 第二个包围盒 + * @return 最短距离 + * + * @ingroup core + */ double distance(const AABB& a, const AABB& b); + +/** + * @brief 两个三角形之间的最短距离 + * + * 处理分离、边相交、面重叠等情况。 + * + * @param a 第一个三角形 + * @param b 第二个三角形 + * @return 最短距离（若相交则为 0） + * + * @ingroup core + */ double distance(const Triangle& a, const Triangle& b); // ── Closest-point queries ──────────────────── + +/** + * @brief 三角形上距给定点最近的点 + * + * 使用重心坐标将投影点截断到三角形内部。 + * 若投影在三角形外，返回最近边或顶点上的点。 + * + * @param p 查询点 + * @param tri 三角形 + * @return 三角形上距离 p 最近的点 + * + * @ingroup core + */ Point3D closest_point(const Point3D& p, const Triangle& tri); + +/** + * @brief 线段上距给定点最近的点 + * + * 通过参数 t = clamp(投影, 0, 1) 实现。 + * + * @param p 查询点 + * @param seg 线段 + * @return 线段上距离 p 最近的点（可能为端点） + * + * @ingroup core + */ Point3D closest_point(const Point3D& p, const Segment3D& seg); + +/** + * @brief 包围盒上距给定点最近的点 + * + * 对每个坐标分量做 clamp 到 [min, max]。 + * + * @param p 查询点 + * @param box 包围盒 + * @return 包围盒边界上（或内部）距离 p 最近的点 + * + * @ingroup core + */ Point3D closest_point(const Point3D& p, const AABB& box); + +/** + * @brief 平面上距给定点最近的点（正交投影） + * + * @param p 查询点 + * @param plane 平面 + * @return 点 p 在平面上的正交投影 + * + * @see Plane::project + * @ingroup core + */ Point3D closest_point(const Point3D& p, const Plane& plane); } // namespace vde::core diff --git a/include/vde/core/icp.h b/include/vde/core/icp.h index 499c3cb..0347af9 100644 --- a/include/vde/core/icp.h +++ b/include/vde/core/icp.h @@ -9,18 +9,62 @@ using foundation::Point3D; using foundation::Vector2D; using foundation::Vector3D; +/** + * @brief ICP（迭代最近点）配准结果 + * + * 包含刚体变换、误差度量、迭代信息和收敛状态。 + * + * @ingroup core + */ struct ICPResult { - Transform3D transform; // Rigid transform aligning source to target - double rms_error; // Root-mean-square error - int iterations; // Number of iterations used + /// 将 source 对齐到 target 的刚体变换 + Transform3D transform; + + /** + * @brief 配准后的均方根误差 + * + * RMS = sqrt(Σ ||target_i - T(source_i)||² / N) + */ + double rms_error; + + /// 实际使用的迭代次数 + int iterations; + + /// 是否在达到最大迭代次数前收敛 bool converged; }; -/// Iterative Closest Point (ICP) for rigid point cloud registration -/// @param source Source point cloud (moved) -/// @param target Target point cloud (fixed) -/// @param max_iter Maximum iterations -/// @param tolerance Convergence tolerance on RMS change +/** + * @brief 迭代最近点（ICP）刚体点云配准 + * + * 使用标准 ICP 算法将 source 点云配准到 target 点云。 + * 每步迭代： + * 1. 对每个 source 点找 target 中的最近邻 + * 2. 用 SVD 求解最优刚体变换 + * 3. 应用变换，检查收敛 + * + * @param source 源点云（将被移动） + * @param target 目标点云（固定参考） + * @param max_iter 最大迭代次数（默认 50） + * @param tolerance 收敛容差：相邻两次 RMS 变化的阈值（默认 1e-6） + * @return 配准结果，包含最终变换和误差 + * + * @note source 和 target 点数可以不同 + * @note 初始对齐越好，收敛越快；若初始偏差太大，可能收敛到局部最优 + * + * @code{.cpp} + * std::vector src = /* ... */; + * std::vector tgt = /* ... */; + * auto result = icp_register(src, tgt, 100, 1e-8); + * if (result.converged) { + * std::cout << "RMS error: " << result.rms_error << "\n"; + * // 应用变换 + * for (auto& p : src) p = result.transform * p; + * } + * @endcode + * + * @ingroup core + */ [[nodiscard]] ICPResult icp_register(const std::vector& source, const std::vector& target, int max_iter = 50, diff --git a/include/vde/core/line.h b/include/vde/core/line.h index b61cd6f..4feb83a 100644 --- a/include/vde/core/line.h +++ b/include/vde/core/line.h @@ -3,49 +3,141 @@ namespace vde::core { +/** + * @brief 三维无限直线 + * + * 由原点位置和单位方向向量表示。 + * 参数方程：P(t) = origin + direction * t。 + * + * @tparam T 标量类型 + * + * @ingroup core + */ template class Line3D { public: + /// 默认构造：原点在 (0,0,0)，方向沿 X 轴 Line3D() = default; + + /** + * @brief 从原点和方向构造直线 + * + * @param origin 直线上一点 + * @param direction 方向向量（自动归一化） + */ Line3D(const Point3D& origin, const Vector3D& direction) : origin_(origin), direction_(direction.normalized()) {} + + /// 获取直线上的一点（原点） [[nodiscard]] const Point3D& origin() const { return origin_; } + /// 获取单位方向向量 [[nodiscard]] const Vector3D& direction() const { return direction_; } + + /** + * @brief 参数方程求点 + * + * @param t 参数值 + * @return origin + direction * t + */ [[nodiscard]] Point3D point_at(T t) const { return origin_ + direction_ * t; } + private: Point3D origin_{0,0,0}; Vector3D direction_{1,0,0}; }; +/** + * @brief 三维线段 + * + * 由两个端点定义的有界直线段。 + * + * @tparam T 标量类型 + * + * @ingroup core + */ template class Segment3D { public: + /// 默认构造：两个端点均在原点 Segment3D() = default; + + /** + * @brief 从两个端点构造线段 + * + * @param p0 起点 + * @param p1 终点 + */ Segment3D(const Point3D& p0, const Point3D& p1) : p0_(p0), p1_(p1) {} + + /// 获取起点 [[nodiscard]] const Point3D& p0() const { return p0_; } + /// 获取终点 [[nodiscard]] const Point3D& p1() const { return p1_; } + + /** + * @brief 方向向量（未归一化） + * @return p1 - p0 + */ [[nodiscard]] Vector3D direction() const { return p1_ - p0_; } + + /** + * @brief 线段长度 + * @return ||p1 - p0|| + */ [[nodiscard]] T length() const { return direction().norm(); } + private: Point3D p0_{0,0,0}, p1_{0,0,0}; }; +/** + * @brief 三维射线 + * + * 由原点和方向定义的半无限直线（t ≥ 0）。 + * 参数方程：P(t) = origin + direction * t, t ≥ 0。 + * + * @tparam T 标量类型 + * + * @ingroup core + */ template class Ray3D { public: + /// 默认构造：原点在 (0,0,0)，方向沿 X 轴 Ray3D() = default; + + /** + * @brief 从原点和方向构造射线 + * + * @param origin 射线起点 + * @param direction 方向向量（自动归一化） + */ Ray3D(const Point3D& origin, const Vector3D& direction) : origin_(origin), direction_(direction.normalized()) {} + + /// 获取射线起点 [[nodiscard]] const Point3D& origin() const { return origin_; } + /// 获取单位方向向量 [[nodiscard]] const Vector3D& direction() const { return direction_; } + + /** + * @brief 参数方程求点（t ≥ 0） + * + * @param t 参数值（应为非负） + * @return origin + direction * t + */ [[nodiscard]] Point3D point_at(T t) const { return origin_ + direction_ * t; } + private: Point3D origin_{0,0,0}; Vector3D direction_{1,0,0}; }; +/// 双精度三维直线 using Line3Dd = Line3D; +/// 双精度三维线段 using Segment3Dd = Segment3D; +/// 双精度三维射线 using Ray3Dd = Ray3D; } // namespace vde::core diff --git a/include/vde/core/plane.h b/include/vde/core/plane.h index d433f81..d6aaf41 100644 --- a/include/vde/core/plane.h +++ b/include/vde/core/plane.h @@ -3,19 +3,84 @@ namespace vde::core { +/** + * @brief 三维平面 + * + * 用 Hessian 法式表示：normal · p + d = 0。 + * 其中 normal 为单位法向量，d 为原点到平面的有向距离的负值。 + * + * @tparam T 标量类型（通常为 double） + * + * @ingroup core + */ template class Plane { public: + /// 默认构造：法线 (0,0,1)，d=0（XY 平面） Plane() = default; + + /** + * @brief 从点和法线构造平面 + * + * @param point 平面上一点 + * @param normal 法向量（自动归一化） + */ Plane(const Point3D& point, const Vector3D& normal) : normal_(normal.normalized()), d_(-normal_.dot(point)) {} + + /** + * @brief 从三点构造平面 + * + * @param a,b,c 平面上不共线的三点（右手定则决定法线方向） + * + * @pre 三点不共线 + */ Plane(const Point3D& a, const Point3D& b, const Point3D& c) : Plane(a, (b - a).cross(c - a)) {} + /// 获取单位法向量 [[nodiscard]] const Vector3D& normal() const { return normal_; } + + /** + * @brief 获取平面方程常数 d + * + * 平面方程为 normal · x + d = 0。 + * d 是原点到平面有向距离的负值。 + * + * @return d 值 + */ [[nodiscard]] T d() const { return d_; } + + /** + * @brief 有向点到平面距离 + * + * 正表示点在法线方向一侧，负表示相反侧。 + * + * @param p 测试点 + * @return normal · p + d + * + * @code{.cpp} + * Plane plane(Point3D(0,0,0), Vector3D::UnitZ()); + * double sd = plane.signed_distance(Point3D(1,2,3)); + * // sd == 3.0 (点在法线方向，正) + * @endcode + */ [[nodiscard]] T signed_distance(const Point3D& p) const { return normal_.dot(p) + d_; } + + /** + * @brief 绝对点到平面距离 + * + * @param p 测试点 + * @return |signed_distance(p)| + */ [[nodiscard]] T distance(const Point3D& p) const { return std::abs(signed_distance(p)); } + + /** + * @brief 点在平面上的正交投影 + * + * @param p 原点 + * @return 投影点（满足在平面上且连线垂直于平面） + */ [[nodiscard]] Point3D project(const Point3D& p) const { return p - normal_ * signed_distance(p); } @@ -25,6 +90,7 @@ private: T d_{0}; }; +/// 双精度三维平面 using Plane3D = Plane; } // namespace vde::core diff --git a/include/vde/core/point.h b/include/vde/core/point.h index 7f2680d..6fd5ba5 100644 --- a/include/vde/core/point.h +++ b/include/vde/core/point.h @@ -2,20 +2,52 @@ #include "vde/foundation/math_types.h" #include "vde/foundation/tolerance.h" +/** + * @defgroup core 核心几何模块 + * + * 核心模块提供 VDE 引擎的几何操作 API： + * - 基础几何体（AABB、直线、射线、线段、平面、三角形、多边形） + * - 几何计算（距离、最近点、凸包、变换） + * - 点云配准（ICP） + * - Voronoi 图 + * + * @{ + */ + namespace vde::core { +/** + * @brief 通用 N 维点类型（重新导出，便于 core 模块使用） + * + * @tparam T 标量类型 + * @tparam Dim 维度 + */ template using Point = foundation::Point; +/// 双精度二维点（重新导出） using Point2D = foundation::Point2D; +/// 双精度三维点（重新导出） using Point3D = foundation::Point3D; +/// 单精度二维点（重新导出） using Point2f = foundation::Point2f; +/// 单精度三维点（重新导出） using Point3f = foundation::Point3f; +/** + * @brief 通用 N 维向量类型（重新导出） + * + * @tparam T 标量类型 + * @tparam Dim 维度 + */ template using Vector = foundation::Vector; +/// 双精度二维向量（重新导出） using Vector2D = foundation::Vector2D; +/// 双精度三维向量（重新导出） using Vector3D = foundation::Vector3D; } // namespace vde::core + +/** @} */ // end of core group diff --git a/include/vde/core/polygon.h b/include/vde/core/polygon.h index df8843f..4448068 100644 --- a/include/vde/core/polygon.h +++ b/include/vde/core/polygon.h @@ -6,44 +6,153 @@ namespace vde::core { +/** + * @brief 二维平面多边形 + * + * 支持面积计算、点包含测试、简化和三角剖分。 + * 顶点按顺序排列（闭合多边形），不存储首尾重复顶点。 + * + * @note 内部假设简单多边形（无自交），非简单多边形的行为未定义 + * @note 所有几何计算基于（假设的）XY 平面 + * + * @ingroup core + */ class Polygon2D { public: + /// 构造空多边形 Polygon2D() = default; + + /** + * @brief 从顶点列表构造多边形 + * + * @param vertices 顶点序列（逆时针或顺时针均可） + */ explicit Polygon2D(std::vector vertices) : vertices_(std::move(vertices)) {} + /// 获取顶点列表 [[nodiscard]] const std::vector& vertices() const { return vertices_; } + /// 获取顶点数量 [[nodiscard]] size_t size() const { return vertices_.size(); } + /// 多边形是否为空（< 3 个顶点视为退化） [[nodiscard]] bool empty() const { return vertices_.empty(); } - /// Signed area (positive = CCW) + /** + * @brief 有向面积（Shoelace 公式） + * + * 正值表示逆时针（CCW），负值表示顺时针（CW），零表示自交或退化。 + * + * @return 有向面积 + * + * @code{.cpp} + * Polygon2D poly({{0,0}, {1,0}, {0,1}}); + * double sa = poly.signed_area(); // > 0 (CCW) + * @endcode + */ [[nodiscard]] double signed_area() const; - /// Absolute area + /** + * @brief 绝对面积 + * @return |signed_area()| + */ [[nodiscard]] double area() const { return std::abs(signed_area()); } - /// Perimeter (sum of edge lengths) + /** + * @brief 周长（各边长之和） + * @return 总边长 + */ [[nodiscard]] double perimeter() const; - /// Area-weighted centroid + /** + * @brief 面积加权质心 + * + * 使用多边形三角剖分的面积加权平均计算质心。 + * 适用于非自交的任意多边形。 + * + * @return 质心坐标 + */ [[nodiscard]] Point2D centroid() const; - /// Axis-aligned bounding box (2D, stored in AABB3D with z=0) + /** + * @brief 二维轴对齐包围盒 + * + * 返回 AABB（z=0），min/max 的 z 分量均为 0。 + * + * @return 包围盒 + */ [[nodiscard]] AABB bounding_box() const; - /// Point-in-polygon test (ray casting) + /** + * @brief 点包含测试（射线法） + * + * 从测试点向右发射水平射线，计算与多边形边的交点数量。 + * 奇数次交点在内部，偶数次在外部。 + * + * @param p 测试点 + * @return true 若点在多边形内部或边界上 + * + * @note 边界上的点也返回 true + */ [[nodiscard]] bool contains(const Point2D& p) const; - /// Is the polygon counter-clockwise? + /** + * @brief 多边形是否为逆时针方向 + * @return signed_area() > 0 + */ [[nodiscard]] bool is_ccw() const { return signed_area() > 0; } - /// Douglas-Peucker simplification. Returns a new simplified polygon. + /** + * @brief Douglas-Peucker 简化 + * + * 用递归的 Douglas-Peucker 算法减少顶点数量， + * 保证简化后多边形与原多边形的最大偏差不超过 tolerance。 + * + * @param tolerance 允许的最大偏差 + * @return 简化后的新多边形 + * + * @code{.cpp} + * // 将 1000 个顶点的多边形简化为约 100 个 + * auto simple = poly.simplify(0.01); + * @endcode + * + * @see triangulate + */ [[nodiscard]] Polygon2D simplify(double tolerance) const; - /// Ear-clipping triangulation. Returns triangles as triples of vertex indices. - /// Assumes a simple polygon (no self-intersections). CCW ordering required. + /** + * @brief 耳切法三角剖分（Ear Clipping） + * + * 将简单多边形剖分为三角形，返回每个三角形的顶点索引三元组。 + * + * @return 三角形列表，每个元素为 {i, j, k} 表示顶点索引 + * + * @pre 多边形为简单多边形（无自交） + * @pre 顶点为逆时针顺序（CCW） + * + * @note 时间复杂度 O(n²)，适用于顶点数 < 1000 的多边形 + * + * @code{.cpp} + * Polygon2D poly(/* CCW vertices */); + * auto tris = poly.triangulate(); + * for (auto& t : tris) { + * // 顶点: poly.vertices()[t[0]], t[1], t[2] + * } + * @endcode + * + * @see simplify + */ [[nodiscard]] std::vector> triangulate() const; - /// Andrew's monotone chain convex hull of a point set. + /** + * @brief Andrew's Monotone Chain 二维凸包 + * + * 计算点集的凸包，返回逆时针排列的凸包多边形。 + * 时间复杂度 O(n log n)。 + * + * @param points 输入点集 + * @return 凸包多边形（CCW，不含共线中间点） + * + * @see convex_hull_2d + */ static Polygon2D convex_hull_2d(const std::vector& points); private: diff --git a/include/vde/core/transform.h b/include/vde/core/transform.h index 6555eca..3fe8dd0 100644 --- a/include/vde/core/transform.h +++ b/include/vde/core/transform.h @@ -4,23 +4,119 @@ namespace vde::core { +/** + * @brief 三维仿射变换类型 + * + * 基于 Eigen::Transform，支持平移、旋转、缩放和它们的组合。 + * 使用 4×4 齐次矩阵表示。 + * + * @ingroup core + */ using Transform3D = Eigen::Transform; +/** + * @brief 创建平移变换 + * + * @param v 平移向量 + * @return 平移变换矩阵 + * + * @code{.cpp} + * auto T = translate(Vector3D(1, 2, 3)); + * Point3D p2 = T * Point3D(0,0,0); // p2 == (1,2,3) + * @endcode + * + * @ingroup core + */ inline Transform3D translate(const Vector3D& v) { return Transform3D(Eigen::Translation(v)); } + +/** + * @brief 创建平移变换（分量形式） + * + * @param x X 方向位移 + * @param y Y 方向位移 + * @param z Z 方向位移 + * @return 平移变换矩阵 + * + * @ingroup core + */ inline Transform3D translate(double x, double y, double z) { return translate(Vector3D(x, y, z)); } + +/** + * @brief 创建旋转变换（轴-角表示） + * + * @param axis 旋转轴（自动归一化） + * @param angle_rad 旋转角度（弧度） + * @return 旋转变换矩阵 + * + * @code{.cpp} + * // 绕 Z 轴旋转 90° + * auto R = rotate(Vector3D::UnitZ(), M_PI / 2); + * @endcode + * + * @ingroup core + */ inline Transform3D rotate(const Vector3D& axis, double angle_rad) { return Transform3D(Eigen::AngleAxis(angle_rad, axis.normalized())); } + +/** + * @brief 绕 X 轴旋转 + * + * @param rad 旋转角度（弧度） + * @return 旋转变换矩阵 + * + * @ingroup core + */ inline Transform3D rotate_x(double rad) { return rotate(Vector3D::UnitX(), rad); } + +/** + * @brief 绕 Y 轴旋转 + * + * @param rad 旋转角度（弧度） + * @return 旋转变换矩阵 + * + * @ingroup core + */ inline Transform3D rotate_y(double rad) { return rotate(Vector3D::UnitY(), rad); } + +/** + * @brief 绕 Z 轴旋转 + * + * @param rad 旋转角度（弧度） + * @return 旋转变换矩阵 + * + * @ingroup core + */ inline Transform3D rotate_z(double rad) { return rotate(Vector3D::UnitZ(), rad); } + +/** + * @brief 创建非均匀缩放变换 + * + * @param sx X 方向缩放因子 + * @param sy Y 方向缩放因子 + * @param sz Z 方向缩放因子 + * @return 缩放变换矩阵 + * + * @note 非均匀缩放可能导致法线方向不准确；使用场景中应考虑法线变换 = (M^{-1})^T + * + * @ingroup core + */ inline Transform3D scale(double sx, double sy, double sz) { return Transform3D(Eigen::Scaling(sx, sy, sz)); } + +/** + * @brief 创建均匀缩放变换 + * + * @param s 各方向统一的缩放因子 + * @return 缩放变换矩阵 + * + * @ingroup core + */ inline Transform3D scale(double s) { return scale(s, s, s); } } // namespace vde::core diff --git a/include/vde/core/triangle.h b/include/vde/core/triangle.h index dbb4f86..ee70fee 100644 --- a/include/vde/core/triangle.h +++ b/include/vde/core/triangle.h @@ -4,24 +4,89 @@ namespace vde::core { +/** + * @brief 三维三角形 + * + * 由三个顶点定义的平面三角形面片。 + * 支持面积、法线、质心、重心坐标包含测试。 + * + * @tparam T 标量类型（通常为 double） + * + * @ingroup core + */ template class Triangle { public: + /// 默认构造：三个顶点均在原点（退化三角形） Triangle() = default; + + /** + * @brief 从三个顶点构造三角形 + * + * @param v0 第一个顶点 + * @param v1 第二个顶点 + * @param v2 第三个顶点 + */ Triangle(const Point3D& v0, const Point3D& v1, const Point3D& v2) : v_{v0, v1, v2} {} + /** + * @brief 获取第 i 个顶点 + * + * @param i 顶点索引（0, 1, 2） + * @return 顶点引用 + */ [[nodiscard]] const Point3D& v(int i) const { return v_[i]; } + + /// 获取所有顶点数组 [[nodiscard]] const std::array& vertices() const { return v_; } + + /** + * @brief 三角形单位法向量 + * + * 使用 (v1 - v0) × (v2 - v0) 的归一化结果。 + * 依右手定则方向。 + * + * @return 单位法向量 + * + * @note 若三角形退化（面积为零），返回零向量 + */ [[nodiscard]] Vector3D normal() const { return (v_[1] - v_[0]).cross(v_[2] - v_[0]).normalized(); } + + /** + * @brief 三角形面积 + * @return 0.5 * ||(v1 - v0) × (v2 - v0)|| + */ [[nodiscard]] T area() const { return (v_[1] - v_[0]).cross(v_[2] - v_[0]).norm() * 0.5; } + + /** + * @brief 三角形质心（几何中心） + * @return (v0 + v1 + v2) / 3 + */ [[nodiscard]] Point3D centroid() const { return (v_[0] + v_[1] + v_[2]) / 3.0; } + + /** + * @brief 重心坐标包含测试 + * + * 将点 p 正交投影到三角形平面上，计算重心坐标并判断是否在三角形内。 + * + * @param p 测试点（三维空间中的任意点） + * @return true 若投影点在三角形内部或边界上 + * + * @note 该测试等效于三维点在三角形上的包含测试（先投影再测试） + * + * @code{.cpp} + * Triangle3D tri({0,0,0}, {1,0,0}, {0,1,0}); + * bool in1 = tri.contains(Point3D(0.2, 0.2, 0)); // true + * bool in2 = tri.contains(Point3D(0.2, 0.2, 5)); // true (投影后判断) + * @endcode + */ [[nodiscard]] bool contains(const Point3D& p) const { Vector3D v0 = v_[2] - v_[0], v1 = v_[1] - v_[0], v2 = p - v_[0]; T d00 = v0.dot(v0), d01 = v0.dot(v1), d11 = v1.dot(v1); @@ -36,7 +101,9 @@ private: std::array v_{}; }; +/// 双精度三维三角形 using Triangle3D = Triangle; +/// 单精度三维三角形 using Triangle3f = Triangle; } // namespace vde::core diff --git a/include/vde/core/voronoi.h b/include/vde/core/voronoi.h index 41d5480..5a3f25a 100644 --- a/include/vde/core/voronoi.h +++ b/include/vde/core/voronoi.h @@ -9,14 +9,56 @@ using foundation::Point3D; using foundation::Vector2D; using foundation::Vector3D; +/** + * @brief Voronoi 单元格 + * + * 表示二维 Voronoi 图中的一个单元格， + * 包含生成点 site 以及该单元格的边界顶点和边。 + * + * @ingroup core + */ struct VoronoiCell { - Point2D site; // Generator point - std::vector vertices; // Cell boundary vertices - std::vector> edges; // Edge indices into vertices + /// 生成点（站点） + Point2D site; + + /** + * @brief 单元格边界顶点（按逆时针顺序排列） + * + * 对于无界单元格，顶点列表可能不闭合。 + */ + std::vector vertices; + + /** + * @brief 单元格边界边 + * + * 每条边表示为 vertices 中的两个索引。 + */ + std::vector> edges; }; -/// 2D Voronoi diagram from Delaunay dual graph -/// Returns cells for each input point +/** + * @brief 二维 Voronoi 图生成（通过 Delaunay 对偶图） + * + * 从输入点集的 Delaunay 三角剖分构造其对偶 Voronoi 图。 + * 每个输入点对应一个 VoronoiCell。 + * + * @param points 输入点集（站点） + * @return 每个站点对应的 VoronoiCell 列表（顺序与输入一致） + * + * @note 对于凸包边界上的站点，对应的 Voronoi 单元格无界，顶点不完全闭合 + * @note 内部使用 Delaunay 三角剖分作为中间步骤 + * + * @code{.cpp} + * std::vector sites = {{0,0}, {1,0}, {0,1}, {1,1}}; + * auto cells = voronoi_2d(sites); + * for (auto& cell : cells) { + * std::cout << "Site: " << cell.site.transpose() << "\n"; + * std::cout << " Vertices: " << cell.vertices.size() << "\n"; + * } + * @endcode + * + * @ingroup core + */ std::vector voronoi_2d(const std::vector& points); } // namespace vde::core diff --git a/include/vde/curves/bezier_curve.h b/include/vde/curves/bezier_curve.h index f28fae5..3ae0807 100644 --- a/include/vde/curves/bezier_curve.h +++ b/include/vde/curves/bezier_curve.h @@ -6,24 +6,88 @@ namespace vde::curves { using core::Point3D; using core::Vector3D; +/** + * @brief 任意阶 Bézier 曲线 + * + * 通过 de Casteljau 递推算法定义的参数多项式曲线。 + * 控制点 cp_ 的数量 = p+1（p 为次数）。 + * 参数域 t ∈ [0, 1]，端点插值首末控制点。 + * + * @ingroup curves + */ class BezierCurve { public: + /** + * @brief 构造 Bézier 曲线 + * @param control_points 控制点序列，点数 = degree + 1，最少 2 个 + * @note 控制点数量决定了曲线阶次；点数越多，曲线越灵活但计算越昂贵 + * @code{.cpp} + * BezierCurve curve({p0, p1, p2, p3}); // 三次 Bézier + * @endcode + */ explicit BezierCurve(std::vector control_points); + /** + * @brief 在参数 t 处求值曲线点 + * @param t 参数值，范围 [0, 1] + * @return 曲线上对应 t 的点坐标 + * @note 使用 de Casteljau 递推算法，O(n²) 复杂度 + */ [[nodiscard]] Point3D evaluate(double t) const; + + /** + * @brief 在参数 t 处求导数 + * @param t 参数值，范围 [0, 1] + * @param order 导数阶数，1 = 一阶（切向量），2 = 二阶（曲率相关） + * @return 导数向量 + * @note 高阶导数通过降阶 Bézier 曲线计算： + * P^{(k)}(t) = n(n-1)...(n-k+1) Σ Δ^k P_i · B_i^{n-k}(t) + */ [[nodiscard]] Vector3D derivative(double t, int order = 1) const; + + /** + * @brief 曲线阶次 + * @return 阶次 = 控制点数 - 1 + */ [[nodiscard]] int degree() const { return static_cast(cp_.size()) - 1; } + + /** + * @brief 参数定义域 + * @return 固定为 [0, 1] + */ [[nodiscard]] std::pair domain() const { return {0.0, 1.0}; } + + /** + * @brief 控制点访问 + * @return 控制点向量的只读引用 + */ [[nodiscard]] const std::vector& control_points() const { return cp_; } - /// Split at parameter t into two curves + /** + * @brief 在 t 处将曲线一分为二 + * @param t 分割参数，范围 [0, 1] + * @return 左段 + 右段两条 Bézier 曲线 + * @note 使用 de Casteljau 分割算法，分割点为三角阵的对角线 + * @code{.cpp} + * auto [left, right] = curve.split(0.5); + * // left: 控制点 = [P0, b1, b2, P(t)] + * // right: 控制点 = [P(t), b3, b4, P3] + * @endcode + */ [[nodiscard]] std::pair split(double t) const; - /// Degree elevation + /** + * @brief 升阶操作 + * @param times 升阶次数，≥ 1 + * @return 升阶后的 Bézier 曲线（控制点数增加 times，形状不变） + * @note 通过反复应用公式 P_i' = i/(n+1) P_{i-1} + (1 - i/(n+1)) P_i + * 升阶不仅增加自由度，也用于不同阶曲线间的兼容操作 + * @see split + */ [[nodiscard]] BezierCurve degree_elevate(int times = 1) const; private: - std::vector cp_; + std::vector cp_; ///< 控制点序列 }; } // namespace vde::curves diff --git a/include/vde/curves/bezier_surface.h b/include/vde/curves/bezier_surface.h index 79f17b1..c04400b 100644 --- a/include/vde/curves/bezier_surface.h +++ b/include/vde/curves/bezier_surface.h @@ -6,20 +6,85 @@ namespace vde::curves { using core::Point3D; using core::Vector3D; +/** + * @brief Bézier 张量积曲面 + * + * 由 (m+1)×(n+1) 控制点网格定义，m = degree_u，n = degree_v。 + * 参数域 (u,v) ∈ [0,1]×[0,1]，四个角点插值网格四角。 + * + * 求值公式：S(u,v) = Σ_i Σ_j P_{ij} · B_i^m(u) · B_j^n(v) + * + * @ingroup curves + */ class BezierSurface { public: + /** + * @brief 构造 Bézier 曲面 + * @param control_grid 控制点网格，grid[i][j] 对应 u 方向第 i 个、v 方向第 j 个控制点 + * @note grid 必须矩整（所有行等长）；grid.size() = degree_u+1, grid[0].size() = degree_v+1 + * @code{.cpp} + * // 3x3 控制网格 → 双二次 Bézier 曲面 + * BezierSurface surf({{p00,p01,p02}, {p10,p11,p12}, {p20,p21,p22}}); + * @endcode + */ BezierSurface(std::vector> control_grid); + /** + * @brief 求曲面点 S(u,v) + * @param u u 参数，范围 [0, 1] + * @param v v 参数，范围 [0, 1] + * @return 曲面上对应 (u,v) 的 3D 坐标 + * @note 张量积求值：先按 v 方向对每行做 de Casteljau，再沿 u 方向对结果做一次 + */ [[nodiscard]] Point3D evaluate(double u, double v) const; + + /** + * @brief u 方向偏导数 ∂S/∂u + * @param u u 参数 + * @param v v 参数 + * @return u 方向切向量 + * @note 通过降阶 Bézier 曲线计算：对 v 求值后沿 u 求导 + */ [[nodiscard]] Vector3D derivative_u(double u, double v) const; + + /** + * @brief v 方向偏导数 ∂S/∂v + * @param u u 参数 + * @param v v 参数 + * @return v 方向切向量 + */ [[nodiscard]] Vector3D derivative_v(double u, double v) const; + + /** + * @brief 单位法向量 N(u,v) = (∂S/∂u × ∂S/∂v) / |…| + * @param u u 参数 + * @param v v 参数 + * @return 归一化法向量（指向曲面的特选侧） + * @note 当两个偏导数平行时返回零向量（奇异点） + */ [[nodiscard]] Vector3D normal(double u, double v) const; + /** + * @brief u 向阶次 + * @return u 方向阶次（控制点行数 - 1） + */ [[nodiscard]] int degree_u() const { return static_cast(cp_.size()) - 1; } + + /** + * @brief v 向阶次 + * @return v 方向阶次（控制点列数 - 1） + */ [[nodiscard]] int degree_v() const { return static_cast(cp_[0].size()) - 1; } private: - std::vector> cp_; + std::vector> cp_; ///< 控制点网格 [u_idx][v_idx] + + /** + * @brief 一维 de Casteljau 求值（内部辅助） + * @param t 参数 + * @param pts 一维控制点序列 + * @return 求值结果 + */ [[nodiscard]] Point3D de_casteljau(double t, const std::vector& pts) const; }; diff --git a/include/vde/curves/bspline_curve.h b/include/vde/curves/bspline_curve.h index 25e7ce6..ae278de 100644 --- a/include/vde/curves/bspline_curve.h +++ b/include/vde/curves/bspline_curve.h @@ -6,32 +6,117 @@ namespace vde::curves { using core::Point3D; using core::Vector3D; +/** + * @brief B 样条曲线（非有理、均匀/非均匀节点） + * + * 由 n+1 个控制点、m+1 个节点向量和阶次 p 定义。 + * 节点向量必须满足 m = n + p + 1。 + * 参数域为 [knots[p], knots[n+1]]，由前 p 个和后 p 个节点修剪。 + * + * 与 Bézier 的主要区别： + * - 局部控制：修改单个控制点仅影响 [t_i, t_{i+p+1}] 区间 + * - 节点插入可引入更多控制点而不改变形状 + * - 通过多重节点可达 C^{p-k} 连续性 + * + * @ingroup curves + */ class BSplineCurve { public: + /** + * @brief 构造 B 样条曲线 + * @param control_points n+1 个控制点 + * @param knots 节点向量，必须为非递减序列 + * @param degree 阶次 p + * @note 节点数量 = 控制点数 + 阶次 + 1；clamped 端点多重度为 p+1 + * @code{.cpp} + * // 三次 B 样条，7 个控制点 + * BSplineCurve curve(cps, {0,0,0,0, 0.25,0.5,0.75, 1,1,1,1}, 3); + * @endcode + */ BSplineCurve(std::vector control_points, std::vector knots, int degree); + /** + * @brief 在参数 t 处求值曲线点 + * @param t 参数值，必须在 domain() 范围内 + * @return 曲线上对应 t 的点坐标 + * @note 使用 de Boor（Cox-de Boor）递推算法 + */ [[nodiscard]] Point3D evaluate(double t) const; + + /** + * @brief 在参数 t 处求导数 + * @param t 参数值 + * @param order 导数阶数，1 = 一切向量 + * @return 导数向量 + * @note B 样条导数公式： + * P^{(k)}(t) = Σ P_i^{(k)} N_{i,p-k}(t)，其中 + * P_i^{(0)} = P_i， + * P_i^{(k)} = p·(P_i^{(k-1)} - P_{i-1}^{(k-1)}) / (t_{i+p-k+1} - t_i) + */ [[nodiscard]] Vector3D derivative(double t, int order = 1) const; + + /** + * @brief 曲线阶次 + * @return 阶次 p + */ [[nodiscard]] int degree() const { return degree_; } + + /** + * @brief 参数定义域 + * @return 有效参数范围 [knots[p], knots[n+1]] + */ [[nodiscard]] std::pair domain() const; + + /** + * @brief 控制点访问 + * @return 控制点数组只读引用 + */ [[nodiscard]] const std::vector& control_points() const { return cp_; } + + /** + * @brief 节点向量访问 + * @return 节点向量只读引用 + */ [[nodiscard]] const std::vector& knots() const { return knots_; } - /// Find knot span index for t + /** + * @brief 查找 t 所在的节点区间 + * @param t 参数值 + * @return 节点区间索引 i 满足 knots[i] ≤ t < knots[i+1] + * @note 使用二分查找，O(log n)；求值和基函数的前置步骤 + * @code{.cpp} + * int span = curve.find_span(0.5); + * auto N = curve.basis_functions(0.5, span); + * @endcode + */ [[nodiscard]] int find_span(double t) const; - /// Evaluate basis functions at t + /** + * @brief 在 t 处求非零基函数值 + * @param t 参数值 + * @param span 节点区间索引，-1 表示自动查找 + * @return N_{span-p}(t), ..., N_{span}(t) 共 p+1 个值 + * @note 仅返回 p+1 个非零基函数；用于端点求值及插值 + * @see find_span + */ [[nodiscard]] std::vector basis_functions(double t, int span = -1) const; - /// Evaluate first derivatives of non-zero basis functions at t - /// Returns dN_{span-degree+i, degree}/du for i = 0..degree + /** + * @brief 在 t 处求非零基函数的一阶导数 + * @param t 参数值 + * @param span 节点区间索引，-1 表示自动查找 + * @return dN_{span-p}/du, ..., dN_{span}/du 共 p+1 个值 + * @note B 样条基函数导数递推： + * N'_{i,p} = p/(t_{i+p}-t_i)·N_{i,p-1} - p/(t_{i+p+1}-t_{i+1})·N_{i+1,p-1} + * @see basis_functions + */ [[nodiscard]] std::vector basis_derivatives(double t, int span = -1) const; private: - std::vector cp_; - std::vector knots_; - int degree_; + std::vector cp_; ///< 控制点序列 + std::vector knots_; ///< 节点向量（非递减） + int degree_; ///< 曲线阶次 }; } // namespace vde::curves diff --git a/include/vde/curves/bspline_surface.h b/include/vde/curves/bspline_surface.h index 2b4952a..e831341 100644 --- a/include/vde/curves/bspline_surface.h +++ b/include/vde/curves/bspline_surface.h @@ -7,26 +7,99 @@ namespace vde::curves { using core::Point3D; using core::Vector3D; +/** + * @brief B 样条张量积曲面（非有理） + * + * 由 (nu+1)×(nv+1) 控制点网格、u/v 向各一个节点向量以及两个阶次定义。 + * 参数域 (u,v) ∈ [knots_u[pu], knots_u[nu+1]] × [knots_v[pv], knots_v[nv+1]]。 + * + * 求值公式：S(u,v) = Σ_i Σ_j P_{ij} · N_{i,pu}(u) · N_{j,pv}(v) + * + * 具有局部控制性：修改单个控制点仅影响相邻 pu×pv 个节点区间的曲面区域。 + * + * @ingroup curves + */ class BSplineSurface { public: + /** + * @brief 构造 B 样条曲面 + * @param control_grid (nu+1)×(nv+1) 控制点网格 + * @param knots_u u 向节点向量，长度 = nu + pu + 1 + * @param knots_v v 向节点向量，长度 = nv + pv + 1 + * @param degree_u u 向阶次 pu + * @param degree_v v 向阶次 pv + * @code{.cpp} + * BSplineSurface surf(grid, + * {0,0,0,0, 0.5, 1,1,1,1}, // u knots (三次，5 个控制点) + * {0,0,0, 0.5, 1,1,1}, // v knots (二次，4 个控制点) + * 3, 2); + * @endcode + */ BSplineSurface(std::vector> control_grid, std::vector knots_u, std::vector knots_v, int degree_u, int degree_v); + /** + * @brief 求曲面点 S(u,v) + * @param u u 参数 + * @param v v 参数 + * @return 曲面上的 3D 点 + * @note 张量积分解：先行后列（或先列后行）做两次 B 样条求值 + */ [[nodiscard]] Point3D evaluate(double u, double v) const; + + /** + * @brief u 向偏导数 ∂S/∂u + * @param u u 参数 + * @param v v 参数 + * @return u 向切向量 + */ [[nodiscard]] Vector3D derivative_u(double u, double v) const; + + /** + * @brief v 向偏导数 ∂S/∂v + * @param u u 参数 + * @param v v 参数 + * @return v 向切向量 + */ [[nodiscard]] Vector3D derivative_v(double u, double v) const; + + /** + * @brief 单位法向量 + * @param u u 参数 + * @param v v 参数 + * @return N(u,v) = (∂S/∂u × ∂S/∂v) / |…| + */ [[nodiscard]] Vector3D normal(double u, double v) const; + /** + * @brief u 向阶次 + * @return pu + */ [[nodiscard]] int degree_u() const { return degree_u_; } + + /** + * @brief v 向阶次 + * @return pv + */ [[nodiscard]] int degree_v() const { return degree_v_; } + + /** + * @brief u 向节点向量 + * @return 只读引用 + */ [[nodiscard]] const std::vector& knots_u() const { return knots_u_; } + + /** + * @brief v 向节点向量 + * @return 只读引用 + */ [[nodiscard]] const std::vector& knots_v() const { return knots_v_; } private: - std::vector> cp_; - std::vector knots_u_, knots_v_; - int degree_u_, degree_v_; + std::vector> cp_; ///< 控制点网格 + std::vector knots_u_, knots_v_; ///< u/v 向节点向量 + int degree_u_, degree_v_; ///< u/v 向阶次 }; } // namespace vde::curves diff --git a/include/vde/curves/nurbs_curve.h b/include/vde/curves/nurbs_curve.h index d8e418d..bb65601 100644 --- a/include/vde/curves/nurbs_curve.h +++ b/include/vde/curves/nurbs_curve.h @@ -7,24 +7,94 @@ namespace vde::curves { using core::Point3D; using core::Vector3D; +/** + * @brief NURBS 曲线（Non-Uniform Rational B-Spline） + * + * B 样条的有理推广，引入权重 w_i 使每个控制点的影响可调。 + * 通过投影变换将齐次坐标中的 B 样条映射到 3D 空间。 + * + * 求值公式：C(t) = Σ N_{i,p}(t) · w_i · P_i / Σ N_{i,p}(t) · w_i + * + * 关键性质： + * - 可精确表示圆锥曲线（圆、椭圆）——Bézier/B 样条无法做到 + * - 权重 w_i 越大，曲线越靠近控制点 P_i + * - 所有权重 = 1 时退化为 B 样条 + * - 所有权重 > 0 时在凸包内 + * + * @ingroup curves + */ class NurbsCurve { public: + /** + * @brief 构造 NURBS 曲线 + * @param control_points 控制点序列 + * @param knots 节点向量 + * @param weights 权重序列，必须与 control_points 等长，各分量 > 0 + * @param degree 阶次 p + * @note 所有权重必须 > 0 以保证分母不为零（非退化） + * @code{.cpp} + * // 用 NURBS 表示四分之一圆（7 个控制点，三次） + * NurbsCurve arc(cps, {0,0,0,0, 0.5, 1,1,1,1}, + * {1, sqrt2/2, 1, sqrt2/2, 1, sqrt2/2, 1}, 3); + * @endcode + */ NurbsCurve(std::vector control_points, std::vector knots, std::vector weights, int degree); + /** + * @brief 在参数 t 处求值 + * @param t 参数值 + * @return 曲线上对应 t 的点坐标 + * @note 内部将 (w·P, w) 提升为齐次坐标的 B 样条，求值后投影回 3D + */ [[nodiscard]] Point3D evaluate(double t) const; + + /** + * @brief 在参数 t 处求导数 + * @param t 参数值 + * @param order 导数阶数 + * @return 导数向量 + * @note 有理导数公式： + * C^{(k)} = (A^{(k)} - Σ_{j=1}^{k} C(k,j) · w^{(j)} · C^{(k-j)}) / w + * 其中 A(t) = Σ N_i(t)·w_i·P_i，w(t) = Σ N_i(t)·w_i + */ [[nodiscard]] Vector3D derivative(double t, int order = 1) const; + + /** + * @brief 曲线阶次 + * @return p + */ [[nodiscard]] int degree() const { return degree_; } + + /** + * @brief 参数定义域 + * @return [knots[degree], knots[n+1]] + */ [[nodiscard]] std::pair domain() const; + + /** + * @brief 控制点访问 + * @return 只读引用 + */ [[nodiscard]] const std::vector& control_points() const { return cp_; } + + /** + * @brief 权重访问 + * @return 只读引用 + */ [[nodiscard]] const std::vector& weights() const { return weights_; } + + /** + * @brief 节点向量访问 + * @return 只读引用 + */ [[nodiscard]] const std::vector& knots() const { return knots_; } private: - std::vector cp_; - std::vector knots_; - std::vector weights_; - int degree_; + std::vector cp_; ///< 控制点 + std::vector knots_; ///< 节点向量 + std::vector weights_; ///< 权重，w_i > 0 + int degree_; ///< 阶次 }; } // namespace vde::curves diff --git a/include/vde/curves/nurbs_operations.h b/include/vde/curves/nurbs_operations.h index 6a9d0fe..2778db2 100644 --- a/include/vde/curves/nurbs_operations.h +++ b/include/vde/curves/nurbs_operations.h @@ -9,79 +9,162 @@ namespace vde::curves { using core::Point3D; using core::Vector3D; -/// Offset a NURBS surface by a constant distance along its normal. -/// Computes the normal at each control point's Greville abscissa and -/// offsets the control point along that normal (Piegl & Tiller approximation). -/// @param surf Input NURBS surface -/// @param distance Offset distance (positive = outward along normal) -/// @return Offset NURBS surface (approximation) +/** + * @brief NURBS 曲面等距偏移 + * + * 基于 Piegl & Tiller 近似方法：在每个控制点的 Greville 横坐标处 + * 计算曲面法向，将控制点沿法向平移 distance。 + * + * @param surf 输入 NURBS 曲面 + * @param distance 偏移距离（正值 = 沿法向外，负值 = 向内） + * @return 偏移后的 NURBS 曲面（近似，阶次与原曲面相同） + * + * @note 这是几何近似，非精确偏移。对于大变曲率曲面，增大控制点数量 + * 可提高精度。精确偏移需要更高阶曲面或细分方法。 + * @note 偏移可能导致自交（distance 超过曲面局部曲率半径时） + * @see trim_surface blend_surfaces + * @ingroup curves + */ [[nodiscard]] NurbsSurface offset_surface(const NurbsSurface& surf, double distance); -/// Trim a NURBS surface to a parameter sub-domain via reparameterization. -/// Maps the old domain [u_min,u_max]×[v_min,v_max] to [0,1]×[0,1] by -/// shifting and scaling knot vectors. Control points and weights are unchanged; -/// the surface geometry is preserved, only the parameter range changes. -/// @param surf Input surface -/// @param u_min, u_max Trimmed u range (must be within original domain) -/// @param v_min, v_max Trimmed v range (must be within original domain) -/// @return Trimmed NURBS surface +/** + * @brief NURBS 曲面参数域裁剪 + * + * 将原始定义域 [u_min,u_max]×[v_min,v_max] 通过节点重参数化 + * 映射到 [0,1]×[0,1]，不改变控制点和权重——仅缩放平移节点向量， + * 曲面几何保持不变。 + * + * @param surf 输入曲面 + * @param u_min,u_max 裁剪后的 u 参数范围（必须在原 domain 内） + * @param v_min,v_max 裁剪后的 v 参数范围（必须在原 domain 内） + * @return 裁剪后的 NURBS 曲面 + * + * @note 裁剪是"软裁剪"——几何不变，仅改变参数化。如需"硬裁剪" + * （沿等参线切除曲面边缘），需结合曲面细分。 + * @code{.cpp} + * auto trimmed = trim_surface(surf, 0.2, 0.8, 0.1, 0.9); + * // trimmed 的参数域为 [0,1]×[0,1]，但几何仅对应原曲面的子区域 + * @endcode + * @ingroup curves + */ [[nodiscard]] NurbsSurface trim_surface(const NurbsSurface& surf, double u_min, double u_max, double v_min, double v_max); -/// Blend between two surfaces along their boundary edges. -/// Extracts boundary curves, offsets them inward by blend_radius along -/// the surface normal, and creates a ruled (lofted) surface between them. -/// @param surf_a First surface -/// @param surf_b Second surface -/// @param edge_a Edge of surf_a to blend from (0=umin, 1=umax, 2=vmin, 3=vmax) -/// @param edge_b Edge of surf_b to blend to (0=umin, 1=umax, 2=vmin, 3=vmax) -/// @param blend_radius Radius of the blend transition -/// @return Blend surface (ruled surface between offset boundary curves) +/** + * @brief 两 NURBS 曲面沿边界边的过渡面 + * + * 提取两个曲面的指定边界曲线，各沿法向向内偏移 blend_radius 后， + * 在两偏移曲线之间创建直纹面作为过渡。 + * + * @param surf_a 第一个曲面 + * @param surf_b 第二个曲面 + * @param edge_a surf_a 的边界边索引（0=u_min, 1=u_max, 2=v_min, 3=v_max） + * @param edge_b surf_b 的边界边索引（0=u_min, 1=u_max, 2=v_min, 3=v_max） + * @param blend_radius 过渡半径 + * @return 过渡 NURBS 曲面（在偏移边界曲线间的直纹面） + * + * @note 过渡面与原始曲面在边界处为 C⁰ 连续（非切线连续） + * @code{.cpp} + * // surf_a 的 u_max 边过渡到 surf_b 的 u_min 边 + * auto blend = blend_surfaces(surf_a, surf_b, + * 1, 0, // edge indices: umax → umin + * 0.5); // blend radius + * @endcode + * @see ruled_surface + * @ingroup curves + */ [[nodiscard]] NurbsSurface blend_surfaces(const NurbsSurface& surf_a, const NurbsSurface& surf_b, int edge_a, int edge_b, double blend_radius); -/// Create a ruled surface between two NURBS curves. -/// The surface is linear (degree=1) in the ruling direction and inherits -/// the curve degree in the longitudinal direction. -/// @param curve_a First curve (v=0 edge of ruled surface) -/// @param curve_b Second curve (v=1 edge of ruled surface) -/// @return Ruled NURBS surface +/** + * @brief 两 NURBS 曲线间的直纹面 + * + * 在两条曲线之间线性插值生成曲面： + * S(u,v) = (1-v)·C_a(u) + v·C_b(u) + * 纵向保持输入曲线的阶次，直纹方向为线性（degree=1）。 + * + * @param curve_a 第一条曲线（v=0 边界） + * @param curve_b 第二条曲线（v=1 边界） + * @return 直纹 NURBS 曲面 + * + * @note 两条曲线需有兼容的阶次和节点结构（内部执行升阶和节点精化） + * @code{.cpp} + * auto ruled = ruled_surface(profile_bottom, profile_top); + * @endcode + * @see coons_patch blend_surfaces + * @ingroup curves + */ [[nodiscard]] NurbsSurface ruled_surface(const NurbsCurve& curve_a, const NurbsCurve& curve_b); -/// Create a Coons patch from four boundary curves. -/// Bilinear interpolation: S(u,v) = (1-u)*C_v0(v) + u*C_v1(v) -/// + (1-v)*C_u0(u) + v*C_u1(u) -/// - corner correction terms. -/// All four curves should have compatible degree and knot structure. -/// @param curve_u0 Boundary at v=0 (runs along u) -/// @param curve_u1 Boundary at v=1 (runs along u) -/// @param curve_v0 Boundary at u=0 (runs along v) -/// @param curve_v1 Boundary at u=1 (runs along v) -/// @return Coons patch NURBS surface +/** + * @brief 由四条边界曲线创建 Coons 曲面 + * + * 双线性混合： + * S(u,v) = (1-u)·C_{v0}(v) + u·C_{v1}(v) + * + (1-v)·C_{u0}(u) + v·C_{u1}(u) + * - 角点修正项 + * + * 四条曲线须围成闭合环且有兼容的阶次/节点结构。 + * + * @param curve_u0 v=0 处边界（沿 u 方向） + * @param curve_u1 v=1 处边界（沿 u 方向） + * @param curve_v0 u=0 处边界（沿 v 方向） + * @param curve_v1 u=1 处边界（沿 v 方向） + * @return Coons 曲面 + * + * @note 仅保证边界插值，不保证内部形状与设计意图一致 + * @code{.cpp} + * auto patch = coons_patch(bottom, top, left, right); + * @endcode + * @see ruled_surface + * @ingroup curves + */ [[nodiscard]] NurbsSurface coons_patch(const NurbsCurve& curve_u0, const NurbsCurve& curve_u1, const NurbsCurve& curve_v0, const NurbsCurve& curve_v1); -/// Extrude a NURBS curve along a direction vector. -/// Creates a degree(d)×1 NURBS surface where d is the curve degree. -/// First row of control points = curve CPs, last row = offset CPs. -/// @param curve Profile curve -/// @param direction Extrusion direction (will be normalized internally) -/// @param length Extrusion length -/// @return Extruded NURBS surface +/** + * @brief 沿方向向量拉伸 NURBS 曲线为曲面 + * + * 创建阶次 (d)×1 的 NURBS 曲面（d = 曲线阶次）。 + * 第一排控制点 = 曲线控制点，最后一排 = 偏移后的控制点。 + * + * @param curve 截面曲线 + * @param direction 拉伸方向（内部归一化） + * @param length 拉伸长度 + * @return 拉伸 NURBS 曲面 + * + * @code{.cpp} + * auto surface = extrude_curve(profile, {0, 0, 1}, 10.0); + * // 沿 Z 轴拉伸 10 单位 + * @endcode + * @ingroup curves + */ [[nodiscard]] NurbsSurface extrude_curve(const NurbsCurve& curve, const Vector3D& direction, double length); -/// Extract a boundary isoparametric curve from a NURBS surface. -/// @param surf Input surface -/// @param edge Which boundary: 0=umin, 1=umax, 2=vmin, 3=vmax -/// @return NURBS curve along the specified boundary +/** + * @brief 提取 NURBS 曲面的等参边界曲线 + * + * 沿曲面的 u/v 参数极值提取边界 NURBS 曲线。 + * + * @param surf 输入曲面 + * @param edge 边界索引：0=u_min, 1=u_max, 2=v_min, 3=v_max + * @return 指定边界上的 NURBS 曲线 + * + * @note 边界曲线继承原曲面该方向的阶次和节点向量 + * @code{.cpp} + * auto umin_curve = extract_boundary_curve(surf, 0); + * auto vmax_curve = extract_boundary_curve(surf, 3); + * @endcode + * @ingroup curves + */ [[nodiscard]] NurbsCurve extract_boundary_curve(const NurbsSurface& surf, int edge); } // namespace vde::curves diff --git a/include/vde/curves/nurbs_surface.h b/include/vde/curves/nurbs_surface.h index 1757126..b331600 100644 --- a/include/vde/curves/nurbs_surface.h +++ b/include/vde/curves/nurbs_surface.h @@ -9,44 +9,154 @@ namespace vde::curves { using core::Point3D; using core::Vector3D; +/** + * @brief NURBS 张量积曲面 + * + * 由 (nu+1)×(nv+1) 控制点网格、权重网格、u/v 节点向量及阶次定义。 + * 求值公式： + * S(u,v) = Σ_i Σ_j N_{i,pu}(u)·N_{j,pv}(v)·w_{ij}·P_{ij} + * / Σ_i Σ_j N_{i,pu}(u)·N_{j,pv}(v)·w_{ij} + * + * 可精确表示球面、圆柱面、锥面等二次曲面。广泛应用于 CAD 数据交换。 + * + * @ingroup curves + */ class NurbsSurface { public: + /** + * @brief 构造 NURBS 曲面 + * @param control_grid 控制点网格 + * @param knots_u u 向节点向量 + * @param knots_v v 向节点向量 + * @param weights 权重网格，维度与 control_grid 一致，所有值 > 0 + * @param degree_u u 向阶次 + * @param degree_v v 向阶次 + * @code{.cpp} + * NurbsSurface sphere(grid, u_knots, v_knots, weights, 3, 3); + * @endcode + */ NurbsSurface(std::vector> control_grid, std::vector knots_u, std::vector knots_v, std::vector> weights, int degree_u, int degree_v); + /** + * @brief 求曲面点 S(u,v) + * @param u u 参数 + * @param v v 参数 + * @return 曲面上对应 (u,v) 的 3D 点 + */ [[nodiscard]] Point3D evaluate(double u, double v) const; + + /** + * @brief u 向偏导数 + * @param u u 参数 + * @param v v 参数 + * @return ∂S/∂u + */ [[nodiscard]] Vector3D derivative_u(double u, double v) const; + + /** + * @brief v 向偏导数 + * @param u u 参数 + * @param v v 参数 + * @return ∂S/∂v + */ [[nodiscard]] Vector3D derivative_v(double u, double v) const; + + /** + * @brief 单位法向量 + * @param u u 参数 + * @param v v 参数 + * @return N(u,v) = normalize(∂S/∂u × ∂S/∂v) + */ [[nodiscard]] Vector3D normal(double u, double v) const; + /** + * @brief u 向阶次 + * @return pu + */ [[nodiscard]] int degree_u() const { return degree_u_; } + + /** + * @brief v 向阶次 + * @return pv + */ [[nodiscard]] int degree_v() const { return degree_v_; } + + /** + * @brief 控制点数量 + * @return {nu+1, nv+1}，即网格的行数和列数 + */ [[nodiscard]] std::array num_control_points() const { return {static_cast(cp_.size()), static_cast(cp_.empty()?0:cp_[0].size())}; } + + /** + * @brief u 向节点向量 + * @return 只读引用 + */ [[nodiscard]] const std::vector& knots_u() const { return knots_u_; } + + /** + * @brief v 向节点向量 + * @return 只读引用 + */ [[nodiscard]] const std::vector& knots_v() const { return knots_v_; } + + /** + * @brief 控制点网格访问 + * @return 只读引用 + */ [[nodiscard]] const std::vector>& control_points() const { return cp_; } + + /** + * @brief 权重网格访问 + * @return 只读引用 + */ [[nodiscard]] const std::vector>& weights() const { return weights_; } - /// Create a ruled surface between two NURBS curves + /** + * @brief 由两条 NURBS 曲线创建直纹面 + * @param c1 起始曲线（v=0 边界） + * @param c2 终止曲线（v=1 边界） + * @return 直纹 NURBS 曲面 + * @note 两条曲线需具有兼容的阶次和节点结构（内部做升阶和节点精化对齐） + * @see blend_surfaces ruled_surface + */ static NurbsSurface ruled(const NurbsCurve& c1, const NurbsCurve& c2); - /// Create a surface of revolution + /** + * @brief 由截面曲线旋转生成旋转曲面 + * @param profile 截面 NURBS 曲线 + * @param axis_origin 旋转轴上一点 + * @param axis_dir 旋转轴方向 + * @param angle_rad 旋转角度（弧度），≤ 2π + * @return 旋转 NURBS 曲面 + * @note 使用 NURBS 表示圆弧的权重方案；完整旋转生成闭曲面 + * @code{.cpp} + * auto rev = NurbsSurface::revolve(profile, origin, {0,0,1}, M_PI); + * @endcode + */ static NurbsSurface revolve(const NurbsCurve& profile, const Point3D& axis_origin, const Vector3D& axis_dir, double angle_rad); - /// Tessellate to triangle mesh + /** + * @brief 将 NURBS 曲面三角化 + * @param res_u u 方向采样分辨率 + * @param res_v v 方向采样分辨率 + * @return {顶点数组, 三角形索引数组}，每个三角形为 {i0,i1,i2} + * @note 生成 res_u×res_v 均匀网格，每个单元格拆为两个三角形 + * @see tessellate + */ [[nodiscard]] std::pair, std::vector>> tessellate(int res_u, int res_v) const; private: - std::vector> cp_; - std::vector knots_u_, knots_v_; - std::vector> weights_; - int degree_u_, degree_v_; + std::vector> cp_; ///< 控制点网格 + std::vector knots_u_, knots_v_; ///< u/v 向节点向量 + std::vector> weights_; ///< 权重网格 + int degree_u_, degree_v_; ///< u/v 向阶次 }; } // namespace vde::curves diff --git a/include/vde/curves/tessellation.h b/include/vde/curves/tessellation.h index aca5dc2..c298eef 100644 --- a/include/vde/curves/tessellation.h +++ b/include/vde/curves/tessellation.h @@ -8,23 +8,58 @@ namespace vde::curves { using core::Point3D; using core::Vector3D; -/// Tessellate a parametric surface into a uniform triangle mesh -/// @param eval Function(u,v) -> Point3D -/// @param res_u, res_v Number of samples in u/v directions -/// @return Pairs of (vertices, triangle indices) +/** + * @brief 均匀参数化曲面三角化 + * + * 在参数域 [0,1]×[0,1] 上生成 res_u×res_v 均匀采样点， + * 调用 eval(u,v) 求值后将每个四边形拆为两个三角形。 + * + * @param eval 曲面求值函数 eval(u,v) -> Point3D，参数 (u,v) ∈ [0,1]×[0,1] + * @param res_u u 方向子分段数（≥ 1） + * @param res_v v 方向子分段数（≥ 1） + * @return {顶点数组, 三角形索引数组}，三角形为 {i0,i1,i2} + * + * @note 总顶点数 = (res_u+1)×(res_v+1)，三角形数 = 2×res_u×res_v + * @code{.cpp} + * auto [verts, tris] = tessellate( + * [&](double u, double v) { return surface.evaluate(u,v); }, + * 20, 20); + * @endcode + * @see adaptive_tessellate + * @ingroup curves + */ std::pair, std::vector>> tessellate(const std::function& eval, int res_u, int res_v); -/// Adaptive tessellation based on normal-angle deviation -/// Recursively subdivides the parameter domain when the angular deviation -/// between corner normals exceeds the threshold. -/// @param eval Function(u,v) -> Point3D -/// @param normal_fn Function(u,v) -> Vector3D (unit normal) -/// @param min_res Minimum subdivisions per direction (1 = single quad) -/// @param max_depth Maximum recursion depth (limits total subdivisions) -/// @param angle_threshold_deg Maximum normal-angle deviation in degrees -/// @return Pairs of (vertices, triangle indices) +/** + * @brief 自适应曲面三角化 + * + * 基于法向角度偏差的自适应细分策略： + * 从 min_res×min_res 初始网格开始，对每个四边形检查四个角点的法向偏差。 + * 若任意相邻角点间的法向夹角超过 angle_threshold_deg，则将该四边形四分细分。 + * + * 相比均匀 tessellate() 的优点： + * - 曲率大的区域自动增加采样密度 + * - 平坦区域保持粗糙网格，节省面数 + * + * @param eval 曲面求值函数 eval(u,v) -> Point3D + * @param normal_fn 单位法向函数 normal_fn(u,v) -> Vector3D + * @param min_res 初始每向子分段数（≥ 1），默认 2 + * @param max_depth 最大递归深度，限制总面数，默认 6 + * @param angle_threshold_deg 法向角度偏差阈值（度），默认 5.0 + * @return {顶点数组, 三角形索引数组} + * + * @note 最终面数上限约为 2·min_res²·4^max_depth（实际远小于此上限） + * @code{.cpp} + * auto [verts, tris] = adaptive_tessellate( + * [&](double u, double v) { return surf.evaluate(u,v); }, + * [&](double u, double v) { return surf.normal(u,v); }, + * 2, 5, 3.0); // 3° 偏差阈值，最大深度 5 + * @endcode + * @see tessellate + * @ingroup curves + */ std::pair, std::vector>> adaptive_tessellate( const std::function& eval, diff --git a/include/vde/foundation/error_codes.h b/include/vde/foundation/error_codes.h index 6446e62..b7ac364 100644 --- a/include/vde/foundation/error_codes.h +++ b/include/vde/foundation/error_codes.h @@ -3,14 +3,30 @@ namespace vde::foundation { +/** + * @brief 全局错误码枚举 + * + * VDE 引擎中所有函数可能返回的错误状态码。 + * 成功用 kSuccess(0) 表示，非零值表示各类错误。 + * + * @ingroup foundation + */ enum class ErrorCode : int32_t { + /// 操作成功完成，无错误 kSuccess = 0, + /// 传入参数无效（空指针、越界、非法值等） kInvalidArgument, + /// 输入几何数据不合法（自交、非流形、退化等） kInvalidGeometry, + /// 退化情况（共线、共面等导致无法确定唯一解） kDegenerateCase, + /// 功能尚未实现 kNotImplemented, + /// 数值计算误差超出容差范围 kNumericalError, + /// 内存分配失败 kOutOfMemory, + /// 内部逻辑错误（不应出现的情况） kInternalError, }; diff --git a/include/vde/foundation/exact_predicates.h b/include/vde/foundation/exact_predicates.h index fa9f34e..24c3503 100644 --- a/include/vde/foundation/exact_predicates.h +++ b/include/vde/foundation/exact_predicates.h @@ -3,32 +3,147 @@ namespace vde::foundation::exact { +/** + * @brief 二维定向测试结果 + * + * 判断三点构成的三角形的定向（顺时针 / 逆时针 / 共线）。 + * + * @ingroup foundation + */ enum class Orientation { Clockwise, CounterClockwise, Collinear }; + +/** + * @brief 圆测试结果 + * + * 判断点 d 相对于三角形 abc 外接圆的位置。 + * + * @ingroup foundation + */ enum class CircleTest { Inside, On, Outside }; -/// Adaptive-precision 2D orientation test (Shewchuk) +/** + * @brief 自适应精度二维定向测试（Shewchuk 算法） + * + * 判断点 a, b, c 的定向。使用自适应浮点精度： + * 先用双精度计算行列式，若结果不可靠则逐步提高精度，直至确定符号。 + * + * @param a 第一个点 + * @param b 第二个点 + * @param c 第三个点 + * @return Orientation::Clockwise 若 a→b→c 顺时针 + * @return Orientation::CounterClockwise 若 a→b→c 逆时针 + * @return Orientation::Collinear 若三点共线 + * + * @note Delaunay 三角剖分、凸包计算等几何算法的核心原语 + * @warning 输入点坐标必须为有限值，NaN/Inf 会导致未定义行为 + * + * @code{.cpp} + * Point2D a(0,0), b(1,0), c(0,1); + * auto ori = exact::orient_2d(a, b, c); + * // ori == Orientation::CounterClockwise + * @endcode + * + * @see orient_3d, in_circle + * @ingroup foundation + */ Orientation orient_2d(const Point2D& a, const Point2D& b, const Point2D& c); -/// Adaptive-precision 3D orientation test +/** + * @brief 自适应精度三维定向测试 + * + * 判断点 a, b, c, d 的定向（即点 d 在平面 abc 的哪一侧）。 + * 等价于计算四面体 abcd 的有向体积的符号。 + * + * @param a 第一个点 + * @param b 第二个点 + * @param c 第三个点 + * @param d 第四个点（测试点） + * @return Orientation::CounterClockwise 若 d 在平面 abc 上方（按右手法则） + * @return Orientation::Clockwise 若 d 在平面 abc 下方 + * @return Orientation::Collinear 若四点共面 + * + * @note 三维 Delaunay 剖分和凸包计算的核心原语 + * + * @code{.cpp} + * Point3D a(0,0,0), b(1,0,0), c(0,1,0), d(0,0,1); + * auto ori = exact::orient_3d(a, b, c, d); + * // ori == Orientation::CounterClockwise + * @endcode + * + * @see orient_2d, in_circle + * @ingroup foundation + */ Orientation orient_3d(const Point3D& a, const Point3D& b, const Point3D& c, const Point3D& d); -/// Adaptive-precision in-circle test +/** + * @brief 自适应精度圆内测试（In-Circle Test） + * + * 判断点 d 是否在三角形 abc 的外接圆内部、边界上或外部。 + * 等价于计算 4×4 行列式的符号。 + * + * @param a 三角形第一个顶点 + * @param b 三角形第二个顶点 + * @param c 三角形第三个顶点 + * @param d 待测试点 + * @return CircleTest::Inside 若 d 在三角形外接圆内部 + * @return CircleTest::On 若 d 恰好在三角形外接圆上 + * @return CircleTest::Outside 若 d 在三角形外接圆外部 + * + * @note Delaunay 三角剖分的空圆性质检查核心原语 + * @pre 三点 a、b、c 不共线 + * + * @code{.cpp} + * Point2D a(0,0), b(2,0), c(0,2); + * Point2D d(0.5, 0.5); // 在外接圆内部 + * auto result = exact::in_circle(a, b, c, d); + * // result == CircleTest::Inside + * @endcode + * + * @see orient_2d, orient_3d + * @ingroup foundation + */ CircleTest in_circle(const Point2D& a, const Point2D& b, const Point2D& c, const Point2D& d); } // namespace vde::foundation::exact -/// Adaptive double-precision predicates (default backend) +/** + * @brief 自适应双精度谓词后端（默认） + * + * 封装 exact 命名空间下的自适应精度谓词函数， + * 作为 Predicates 的默认后端实现。 + * + * @ingroup foundation + */ namespace vde::foundation { struct AdaptiveDouble { + /** + * @brief 二维定向测试 + * @param a,b,c 三个二维点 + * @return 定向结果 + * @see exact::orient_2d + */ static exact::Orientation orient_2d(const Point2D& a, const Point2D& b, const Point2D& c) { return exact::orient_2d(a, b, c); } + /** + * @brief 三维定向测试 + * @param a,b,c,d 四个三维点 + * @return 定向结果 + * @see exact::orient_3d + */ static exact::Orientation orient_3d(const Point3D& a, const Point3D& b, const Point3D& c, const Point3D& d) { return exact::orient_3d(a, b, c, d); } + /** + * @brief 圆内测试 + * @param a,b,c 三角形三顶点 + * @param d 待测试点 + * @return 圆测试结果 + * @see exact::in_circle + */ static exact::CircleTest in_circle(const Point2D& a, const Point2D& b, const Point2D& c, const Point2D& d) { return exact::in_circle(a, b, c, d); diff --git a/include/vde/foundation/exact_predicates_gmp.h b/include/vde/foundation/exact_predicates_gmp.h index 90c907c..2c5bd98 100644 --- a/include/vde/foundation/exact_predicates_gmp.h +++ b/include/vde/foundation/exact_predicates_gmp.h @@ -6,11 +6,39 @@ namespace vde::foundation::exact { -/// GMP-based exact predicates — zero rounding error, arbitrary precision +/** + * @brief GMP 精确算术谓词后端 + * + * 基于 GNU MP（多精度有理数）的精确几何谓词。 + * 使用高精度有理数运算，零舍入误差，适用于对鲁棒性要求极高的场景。 + * + * @note 仅在定义了 VDE_USE_GMP 宏且链接了 GMP 库时可用 + * @warning GMP 版本比自适应双精度版本慢约 10-50 倍，仅在双精度不足以确定结果时使用 + * + * @ingroup foundation + */ struct GmpExact { + /// GMP 有理数类型 using Rational = mpq_class; - /// orient_2d with exact rational arithmetic + /** + * @brief 精确有理数二维定向测试 + * + * 使用 GMP 有理数精确计算三点定向，零误差。 + * + * @param a 第一个点 + * @param b 第二个点 + * @param c 第三个点 + * @return 定向结果（保证正确，无舍入误差） + * + * @code{.cpp} + * // 对近共线点仍能给出确定结果 + * Point2D a(0, 0), b(1, 1e-15), c(2, 2e-15); + * auto ori = GmpExact::orient_2d(a, b, c); + * @endcode + * + * @see GmpExact::orient_3d, GmpExact::in_circle + */ static Orientation orient_2d(const Point2D& a, const Point2D& b, const Point2D& c) { Rational ax(a.x()), ay(a.y()); Rational bx(b.x()), by(b.y()); @@ -21,7 +49,17 @@ struct GmpExact { s < 0 ? Orientation::Clockwise : Orientation::Collinear; } - /// orient_3d with exact rational arithmetic + /** + * @brief 精确有理数三维定向测试 + * + * 判断点 d 相对于平面 abc 的位置，使用有理数精确运算。 + * + * @param a,b,c 定义平面的三点 + * @param d 待测试点 + * @return 定向结果（零误差） + * + * @see GmpExact::orient_2d + */ static Orientation orient_3d(const Point3D& a, const Point3D& b, const Point3D& c, const Point3D& d) { Rational ax(a.x()), ay(a.y()), az(a.z()); @@ -42,7 +80,17 @@ struct GmpExact { s < 0 ? Orientation::Clockwise : Orientation::Collinear; } - /// in_circle with exact rational arithmetic + /** + * @brief 精确有理数圆内测试 + * + * 判断点 d 是否在三角形 abc 的外接圆内，零误差。 + * + * @param a,b,c 三角形三顶点 + * @param d 待测试点 + * @return 圆测试结果（零误差） + * + * @see GmpExact::orient_2d + */ static CircleTest in_circle(const Point2D& a, const Point2D& b, const Point2D& c, const Point2D& d) { Rational ax(a.x()), ay(a.y()); diff --git a/include/vde/foundation/format_utils.h b/include/vde/foundation/format_utils.h index 74b7e16..77ca9f0 100644 --- a/include/vde/foundation/format_utils.h +++ b/include/vde/foundation/format_utils.h @@ -8,7 +8,26 @@ namespace vde::foundation { -/// Format double with adequate precision for CAD (12 significant digits) +/** + * @brief 以 CAD 级别精度格式化双精度浮点数 + * + * 输出科学计数法字符串，保证 12 位有效数字，满足 CAD 数据交换精度要求。 + * NaN 被格式化为 "0.0"，Inf/ -Inf 被裁剪为安全边界值。 + * + * @param v 要格式化的双精度数 + * @param precision 有效数字位数（默认 12） + * @return 格式化后的字符串，如 "1.234567890123e+02" + * + * @note 总是确保指数部分带符号（如 "e+00" 而非 "e00"），便于解析器统一处理 + * + * @code{.cpp} + * auto s = fmt_double(3.14159265358979); + * // s == "3.141592653590e+00" + * @endcode + * + * @see fmt_iges_real, fmt_point + * @ingroup foundation + */ inline std::string fmt_double(double v, int precision = 12) { if (std::isnan(v)) return "0.0"; if (std::isinf(v)) return v > 0 ? "1e30" : "-1e30"; @@ -25,7 +44,22 @@ inline std::string fmt_double(double v, int precision = 12) { return s; } -/// Format for IGES: uses D exponent +/** + * @brief IGES 格式实数输出 + * + * 与 fmt_double() 相同，但使用 IGES 标准的 'D' 指数标记代替 'E'。 + * + * @param v 要格式化的双精度数 + * @return IGES 格式字符串，如 "3.141592653590D+00" + * + * @code{.cpp} + * auto s = fmt_iges_real(3.14159265358979); + * // s == "3.141592653590D+00" + * @endcode + * + * @see fmt_double + * @ingroup foundation + */ inline std::string fmt_iges_real(double v, int /*precision*/ = 12) { auto s = fmt_double(v, 12); auto pos = s.find('e'); @@ -33,14 +67,41 @@ inline std::string fmt_iges_real(double v, int /*precision*/ = 12) { return s; } -/// Format a Point3D as comma-separated reals +/** + * @brief 格式化三维点为逗号分隔字符串 + * + * 输出 "x,y,z" 形式，常用于 IGES/STEP/OBJ 等文本格式的顶点坐标行。 + * + * @param p 三维点 + * @param use_d 若为 true，使用 IGES 风格 D 指数标记 + * @return 逗号分隔的坐标字符串 + * + * @code{.cpp} + * Point3D p(1.0, 2.5, -3.14); + * auto s = fmt_point(p); + * // s == "1.000000000000e+00,2.500000000000e+00,-3.140000000000e+00" + * @endcode + * + * @see fmt_double, fmt_iges_real + * @ingroup foundation + */ inline std::string fmt_point(const core::Point3D& p, bool use_d = false) { using FmtFn = std::string (*)(double, int); FmtFn fn = use_d ? static_cast(fmt_iges_real) : fmt_double; return fn(p.x(), 12) + "," + fn(p.y(), 12) + "," + fn(p.z(), 12); } -/// Pad a string to exactly 72 characters with spaces +/** + * @brief 将字符串填充到恰好 72 个字符 + * + * IGES 格式中每条数据行固定为 72 字符宽。不足用空格补齐，超出则截断。 + * + * @param s 输入字符串 + * @return 恰好 72 字符的字符串（不足右侧补空格，超出截断） + * + * @see format_iges_line + * @ingroup foundation + */ inline std::string pad72(const std::string& s) { std::string result = s; if (result.size() > 72) result.resize(72); @@ -48,7 +109,24 @@ inline std::string pad72(const std::string& s) { return result; } -/// Format an IGES section line: data + section marker (80 columns total) +/** + * @brief 格式化 IGES 段落行 + * + * 生成一条完整的 IGES 80 列记录行：72 列数据 + 8 列序号（section_char + 7 位数字）。 + * + * @param data 数据内容（最多 72 字符，超出截断） + * @param section_char 段落标识符（如 'P' 参数段, 'D' 目录段） + * @param seq 行序号（1-based） + * @return 完整 IGES 行（含换行符），共 81 字符 + * + * @code{.cpp} + * auto line = format_iges_line("110,1.0,2.0,3.0;", 'P', 1); + * // 输出: "110,1.0,2.0,3.0; P0000001\n" + * @endcode + * + * @see pad72 + * @ingroup foundation + */ inline std::string format_iges_line(const std::string& data, char section_char, int seq) { std::string line = pad72(data); char seq_buf[9]; diff --git a/include/vde/foundation/interval.h b/include/vde/foundation/interval.h index 60f254b..effc54f 100644 --- a/include/vde/foundation/interval.h +++ b/include/vde/foundation/interval.h @@ -5,33 +5,94 @@ namespace vde::foundation { -/// 双精度区间算术：存储 [lo, hi] 边界，用于验证计算结果可靠性 +/** + * @brief 双精度区间算术类 + * + * 存储区间 [lo, hi] 的上下边界，支持区间算术运算和确定性比较。 + * 主要用于验证浮点计算的可靠性：若区间结果不跨越零，则可确定符号； + * 若跨越零，则说明双精度不足以确定结果，需要更高精度或精确算术。 + * + * 区间运算遵循保守原则（结果区间包含所有可能的真实值）， + * 即 out = [min(所有组合), max(所有组合)]。 + * + * @ingroup foundation + */ class Interval { public: double lo, hi; + /// 默认构造：空区间 [0, 0] Interval() : lo(0), hi(0) {} + + /** + * @brief 从单值构造精确区间（退化为点） + * @param v 区间值，[v, v] + */ Interval(double v) : lo(v), hi(v) {} + + /** + * @brief 从上下界构造区间 + * + * 自动交换以确保 lo ≤ hi。 + * + * @param l 下界 + * @param h 上界 + */ Interval(double l, double h) : lo(std::min(l,h)), hi(std::max(l,h)) {} - /// 区间中点 + /** + * @brief 区间中点 + * @return (lo + hi) / 2 + */ [[nodiscard]] double mid() const { return (lo + hi) * 0.5; } - /// 区间宽度 + /** + * @brief 区间宽度 + * @return hi - lo + */ [[nodiscard]] double width() const { return hi - lo; } - /// 区间是否包含值 + /** + * @brief 判断区间是否包含给定值 + * @param v 测试值 + * @return true 若 lo ≤ v ≤ hi + */ [[nodiscard]] bool contains(double v) const { return lo <= v && v <= hi; } - /// 是否精确（宽度 < eps） + /** + * @brief 判断区间是否足够窄（收敛） + * @param eps 宽度阈值（默认 1e-12） + * @return true 若区间宽度 < eps + */ [[nodiscard]] bool is_exact(double eps = 1e-12) const { return width() < eps; } - /// 区间是否与另一个区间重叠 + /** + * @brief 判断两个区间是否有重叠 + * @param o 另一个区间 + * @return true 若存在公共点 + */ [[nodiscard]] bool overlaps(const Interval& o) const { return !(hi < o.lo || o.hi < lo); } - /// 确定性比较 + /** + * @brief 确定性区间比较 + * + * 仅当两个区间无重叠时才能确定大小关系。 + * + * @param o 另一个区间 + * @return -1 若确定 this < o + * @return 1 若确定 this > o + * @return 0 若区间重叠，无法确定 + * + * @code{.cpp} + * Interval a(1.0, 2.0), b(3.0, 4.0); + * int c = a.cmp(b); // c == -1 (确定 a < b) + * + * Interval c(1.0, 3.0), d(2.0, 4.0); + * int r = c.cmp(d); // r == 0 (重叠，无法确定) + * @endcode + */ [[nodiscard]] int cmp(const Interval& o) const { if (hi < o.lo) return -1; // 确定 this < o if (lo > o.hi) return 1; // 确定 this > o @@ -39,16 +100,42 @@ public: } // ── 算术运算 ── + + /** + * @brief 区间加法 + * @return [this.lo + o.lo, this.hi + o.hi] + */ Interval operator+(const Interval& o) const { return {lo + o.lo, hi + o.hi}; } + + /** + * @brief 区间减法 + * @return [this.lo - o.hi, this.hi - o.lo]（保守估计） + */ Interval operator-(const Interval& o) const { return {lo - o.hi, hi - o.lo}; } + + /** + * @brief 区间乘法 + * + * 取四项乘积的最小值和最大值作为结果区间。 + * + * @return [min(a*c,a*d,b*c,b*d), max(a*c,a*d,b*c,b*d)] + */ 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})}; } + + /** + * @brief 区间除法 + * + * 若除区间跨越零，返回极大区间表示不确定性。 + * + * @return 结果区间或 [-∞, +∞]（除数为零区间） + */ Interval operator/(const Interval& o) const { if (o.lo <= 0 && o.hi >= 0) return {-std::numeric_limits::max(), std::numeric_limits::max()}; @@ -56,26 +143,60 @@ public: return {std::min({a,b,c,d}), std::max({a,b,c,d})}; } + /// 取负 Interval operator-() const { return {-hi, -lo}; } // ── 数学函数（保守区间）── + + /** + * @brief 区间平方根 + * + * 下界对 0 取 max 以处理负值（保守处理）。 + * + * @return [sqrt(max(0, lo)), sqrt(max(0, hi))] + */ [[nodiscard]] Interval sqrt() const { return {std::sqrt(std::max(0.0, lo)), std::sqrt(std::max(0.0, hi))}; } - /// 从 double 创建区间（含舍入误差） + /** + * @brief 从双精度值创建区间（含舍入误差带） + * + * @param v 中心值 + * @param error 误差半径（默认 0） + * @return [v - error, v + error] + */ static Interval from_double(double v, double error = 0) { return {v - error, v + error}; } - /// 将 double 表示的坐标转为区间（含 ULP 误差） + /** + * @brief 从坐标值创建区间（含 ULP 误差带） + * + * 误差估计为 |v| * 2.22e-16 + 1e-300，即双精度最小单位。 + * + * @param v 坐标值 + * @return 包含浮点舍入误差的区间 + */ static Interval from_coord(double v) { double eps = std::abs(v) * 2.22e-16 + 1e-300; return {v - eps, v + eps}; } }; -/// 行列式的区间验证 +/** + * @brief 二维定向行列式的区间验证 + * + * 用区间算术计算 orient_2d 行列式，结果区间用于判断双精度是否可靠。 + * + * @param ax,ay 点 a 的坐标 + * @param bx,by 点 b 的坐标 + * @param cx,cy 点 c 的坐标 + * @return 行列式的区间值 + * + * @see orient_2d_verified + * @ingroup foundation + */ 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); @@ -84,8 +205,30 @@ inline Interval orient_2d_interval(double ax, double ay, double bx, double by, d return acx * bcy - acy * bcx; } -/// 使用区间算术验证 orient_2d 结果 -/// 返回: -1(CW), 0(无法确定), 1(CCW) +/** + * @brief 使用区间算术验证 orient_2d 结果 + * + * 若双精度计算结果接近零，调用此函数进行区间验证。 + * 若区间不跨越零，则结果可靠；否则需要升级到精确算术。 + * + * @param ax,ay 点 a 的坐标 + * @param bx,by 点 b 的坐标 + * @param cx,cy 点 c 的坐标 + * @return -1 确认 Clockwise + * @return 0 无法确定（需升级到 GMP / 自适应精度） + * @return 1 确认 CounterClockwise + * + * @code{.cpp} + * int result = orient_2d_verified(0.0, 0.0, 1.0, 1e-16, 2.0, 2e-16); + * if (result == 0) { + * // 双精度不可靠，需要精确算术 + * auto ori = GmpExact::orient_2d(a, b, c); + * } + * @endcode + * + * @see orient_2d_interval, Interval::cmp + * @ingroup foundation + */ 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)); diff --git a/include/vde/foundation/io_gltf.h b/include/vde/foundation/io_gltf.h index c3f6c95..a5a77e4 100644 --- a/include/vde/foundation/io_gltf.h +++ b/include/vde/foundation/io_gltf.h @@ -7,24 +7,76 @@ namespace vde::foundation { -/// GLTF export options +/** + * @brief glTF 2.0 导出选项 + * + * 控制网格以 glTF/GLB 格式导出的参数。 + * + * @ingroup foundation + */ struct GltfOptions { - bool binary = true; // GLB (single file) vs glTF (JSON + .bin) - bool include_normals = true; // Compute and export vertex normals - bool include_colors = false; // Per-vertex colors - std::vector vertex_colors; // RGB per vertex (if include_colors) + /// 是否以二进制 GLB 格式（单文件）输出；false 则输出 JSON + .bin 分离格式 + bool binary = true; + + /// 是否计算和导出顶点法线 + bool include_normals = true; + + /// 是否导出逐顶点颜色 + bool include_colors = false; + + /// 逐顶点 RGB 颜色数据（仅当 include_colors = true 时使用） + std::vector vertex_colors; }; -/// Export mesh to glTF 2.0 / GLB -/// Returns true on success +/** + * @brief 导出网格为 glTF 2.0 / GLB 格式 + * + * 将 HalfedgeMesh 写入 glTF 2.0 规范文件。默认输出单文件 GLB 二进制格式， + * 可通过 GltfOptions 控制输出格式和包含的属性。 + * + * @param filepath 输出文件路径（.glb 或 .gltf） + * @param mesh 要导出的半边网格 + * @param options 导出选项（格式、法线、颜色等） + * @return true 导出成功 + * @return false 导出失败（文件无法创建等） + * + * @note 生成的 glTF 符合 2.0 规范，兼容 Blender、three.js、Unreal Engine 等工具 + * + * @code{.cpp} + * mesh::HalfedgeMesh mesh = /* ... */; + * GltfOptions opts; + * opts.binary = true; + * opts.include_normals = true; + * bool ok = write_gltf("output.glb", mesh, opts); + * @endcode + * + * @see write_brep_gltf, GltfOptions + * @ingroup foundation + */ bool write_gltf(const std::string& filepath, const mesh::HalfedgeMesh& mesh, const GltfOptions& options = GltfOptions{}); -/// Export B-Rep model to GLB by tessellating faces -/// @param filepath Output file (.glb recommended) -/// @param model B-Rep model to export -/// @param tessellation_res Resolution for curved surfaces (segments per edge) -/// @return true on success +/** + * @brief 导出 B-Rep 模型为 GLB 格式（通过面片化） + * + * 将 B-Rep 模型的参数曲面离散化为三角网格后导出为 GLB。 + * 曲面通过 tessellation_res 控制离散精度。 + * + * @param filepath 输出文件路径（建议 .glb） + * @param model 要导出的 B-Rep 模型 + * @param tessellation_res 每条边的离散段数（默认 32），越大越精确 + * @return true 导出成功 + * @return false 导出失败 + * + * @code{.cpp} + * brep::BrepModel model = /* ... */; + * // 高精度曲面导出 + * bool ok = write_brep_gltf("model.glb", model, 64); + * @endcode + * + * @see write_gltf + * @ingroup foundation + */ bool write_brep_gltf(const std::string& filepath, const brep::BrepModel& model, int tessellation_res = 32); diff --git a/include/vde/foundation/io_obj.h b/include/vde/foundation/io_obj.h index 5409d35..cbc883d 100644 --- a/include/vde/foundation/io_obj.h +++ b/include/vde/foundation/io_obj.h @@ -5,31 +5,106 @@ namespace vde::foundation { +/** + * @brief OBJ 文件解析结果数据结构 + * + * 存储从 Wavefront OBJ 文件读取的完整几何数据， + * 包括顶点位置、纹理坐标、法线、面信息和材质分组。 + * + * 内部索引均为 0-based，与文件中的 1-based 索引不同。 + * + * @ingroup foundation + */ struct ObjMeshData { + /// 顶点位置数组 std::vector vertices; + /// 纹理坐标数组（可能为空） std::vector texcoords; + /// 顶点法线数组（可能为空） std::vector normals; - /// Vertex indices per face (0-based internally). - /// Always populated regardless of face format. + /** + * @brief 每个面的顶点索引列表（0-based） + * + * 无论文件中面的格式如何（纯顶点、"v/vt"、"v/vt/vn"）， + * 此字段始终被填充。faces[i] 存储第 i 个面的所有顶点索引。 + */ std::vector> faces; - /// Texcoord indices per face (0-based). Empty if the file had no vt data. - /// Same length as faces; each sub-vector same length as corresponding face. + /** + * @brief 每个面的纹理坐标索引列表（0-based） + * + * 与 faces 数组长度相同；若无纹理数据则为空。 + */ std::vector> face_texcoords; - /// Normal indices per face (0-based). Empty if the file had no vn data. + /** + * @brief 每个面的法线索引列表（0-based） + * + * 与 faces 数组长度相同；若无法线数据则为空。 + */ std::vector> face_normals; - /// Material name per face group (appears before the faces that use it). - /// The i-th entry gives the material active for the i-th face. + /** + * @brief 每个面对应的材质名称 + * + * face_materials[i] 给出第 i 个面使用的材质名。 + * 材质名来源于 usemtl 指令，在定义后所有后续面使用该材质。 + */ std::vector face_materials; - /// Helpers + /// 是否有纹理坐标数据 [[nodiscard]] bool has_texcoords() const { return !texcoords.empty(); } }; +/** + * @brief 读取 Wavefront OBJ 文件 + * + * 解析 ASCII OBJ 文件，支持以下 OBJ 特性： + * - v（顶点）、vt（纹理坐标）、vn（法线） + * - f（面）：支持三角面和任意多边形面 + * - usemtl（材质名称）、g（组名）、o（对象名） + * - # 注释行 + * + * @param filepath OBJ 文件路径 + * @return 解析后的 ObjMeshData + * + * @throws std::runtime_error 若文件无法打开或格式错误 + * + * @note 不解析 MTL 材质文件，仅记录 usemtl 指令中的材质名称 + * @note 面索引在内部转换为 0-based；若文件中使用相对索引（负值），自动转换为正向 + * + * @code{.cpp} + * auto data = read_obj("model.obj"); + * std::cout << "Vertices: " << data.vertices.size() << "\n"; + * std::cout << "Faces: " << data.faces.size() << "\n"; + * @endcode + * + * @see write_obj, ObjMeshData + * @ingroup foundation + */ ObjMeshData read_obj(const std::string& filepath); + +/** + * @brief 写入 Wavefront OBJ 文件 + * + * 将 ObjMeshData 以 ASCII OBJ 格式写出，包含顶点、纹理坐标、 + * 法线和面信息。材质名称通过注释嵌入。 + * + * @param filepath 输出文件路径 + * @param data 要写出的网格数据 + * + * @note 输出的面索引为 1-based（OBJ 标准） + * + * @code{.cpp} + * ObjMeshData data; + * // ... 填充 data ... + * write_obj("output.obj", data); + * @endcode + * + * @see read_obj + * @ingroup foundation + */ void write_obj(const std::string& filepath, const ObjMeshData& data); } // namespace vde::foundation diff --git a/include/vde/foundation/io_ply.h b/include/vde/foundation/io_ply.h index 75a07c8..e6e5062 100644 --- a/include/vde/foundation/io_ply.h +++ b/include/vde/foundation/io_ply.h @@ -5,11 +5,51 @@ namespace vde::foundation { -/// Read PLY file (Stanford Polygon Format) -/// Supports ASCII and binary little-endian, vertex + face elements +/** + * @brief 读取 PLY 文件（Stanford Polygon Format） + * + * 解析 PLY 格式文件，支持 ASCII 和二进制（小端序）两种编码。 + * 支持的 PLY 元素：vertex（位置、法线、颜色）、face（顶点索引列表）。 + * + * @param filepath PLY 文件路径 + * @return 解析后的半边网格 + * + * @throws std::runtime_error 若文件无法打开、格式不支持或数据损坏 + * + * @note 自动检测文件头部 magic number 判断 ASCII / 二进制格式 + * @note 二进制格式仅支持小端序（little-endian） + * @note 返回类型为 HalfedgeMesh，内部自动建立半边拓扑 + * + * @code{.cpp} + * auto mesh = read_ply("scan.ply"); + * std::cout << "Vertices: " << mesh.vertex_count() << "\n"; + * std::cout << "Faces: " << mesh.face_count() << "\n"; + * @endcode + * + * @see write_ply + * @ingroup foundation + */ mesh::HalfedgeMesh read_ply(const std::string& filepath); -/// Write PLY file (ASCII format) +/** + * @brief 写入 PLY 文件（ASCII 格式） + * + * 将半边网格以 ASCII PLY 格式写出，包含顶点位置和面信息。 + * + * @param filepath 输出文件路径 + * @param mesh 要导出的半边网格 + * + * @note 输出的 PLY 为 ASCII 编码，易于人工阅读和编辑 + * @note 若需要二进制输出以提高性能，请使用 BinarySerializer + * + * @code{.cpp} + * mesh::HalfedgeMesh mesh = /* ... */; + * write_ply("output.ply", mesh); + * @endcode + * + * @see read_ply, BinarySerializer + * @ingroup foundation + */ void write_ply(const std::string& filepath, const mesh::HalfedgeMesh& mesh); } // namespace vde::foundation diff --git a/include/vde/foundation/io_stl.h b/include/vde/foundation/io_stl.h index 8bbe61d..e4cdc98 100644 --- a/include/vde/foundation/io_stl.h +++ b/include/vde/foundation/io_stl.h @@ -5,18 +5,85 @@ namespace vde::foundation { +/** + * @brief STL 三角形面片 + * + * 表示一个 STL 三角形，包含法向量和三个顶点。 + * 法向量应为单位向量，方向朝外（右手法则）。 + * + * @ingroup foundation + */ struct StlTriangle { + /// 面法向量（单位向量） Vector3D normal; + /// 三个顶点（右手法则 CCW 顺序） Point3D v0, v1, v2; }; -/// Auto-detect format (binary or ASCII) and read. +/** + * @brief 读取 STL 文件（自动检测格式） + * + * 通过检查文件头部的 80 字节判断二进制（含 triangle count） + * 还是 ASCII（以 "solid" 开头）格式，然后选择对应解析器。 + * + * @param filepath STL 文件路径 + * @return 三角形面片列表 + * + * @throws std::runtime_error 若文件无法打开或格式无效 + * + * @note 自动格式检测：二进制 STL 的 80 字节头后跟 4 字节三角形数量 + * @note ASCII STL 以 "solid name" 开头 + * + * @code{.cpp} + * auto tris = read_stl("part.stl"); + * std::cout << "Triangles: " << tris.size() << "\n"; + * @endcode + * + * @see write_stl, write_stl_ascii + * @ingroup foundation + */ std::vector read_stl(const std::string& filepath); -/// Always write binary STL. +/** + * @brief 写入二进制 STL 文件 + * + * 以紧凑的二进制格式写出 STL 三角形列表。 + * 二进制 STL 为业界标准交换格式，文件体积小、读写快。 + * + * @param filepath 输出文件路径 + * @param tris 三角形面片列表 + * + * @note 二进制 STL 格式：80 字节头 + 4 字节三角形数 + 每三角形 50 字节 + * @note 法向量在写入时自动归一化 + * + * @code{.cpp} + * std::vector tris = /* ... */; + * write_stl("output.stl", tris); // 二进制，体积小 + * @endcode + * + * @see read_stl, write_stl_ascii + * @ingroup foundation + */ void write_stl(const std::string& filepath, const std::vector& tris); -/// Write ASCII STL (human-readable). +/** + * @brief 写入 ASCII STL 文件 + * + * 以人类可读的 ASCII 文本格式写出 STL，便于调试和手动编辑。 + * + * @param filepath 输出文件路径 + * @param tris 三角形面片列表 + * + * @warning ASCII STL 文件体积约为二进制 STL 的 5-10 倍，不适合大规模网格 + * + * @code{.cpp} + * std::vector tris = /* ... */; + * write_stl_ascii("debug.stl", tris); // 可读文本，便于调试 + * @endcode + * + * @see read_stl, write_stl + * @ingroup foundation + */ void write_stl_ascii(const std::string& filepath, const std::vector& tris); } // namespace vde::foundation diff --git a/include/vde/foundation/math_types.h b/include/vde/foundation/math_types.h index 79b7ade..8c72e01 100644 --- a/include/vde/foundation/math_types.h +++ b/include/vde/foundation/math_types.h @@ -5,30 +5,79 @@ namespace vde::foundation { +/** + * @brief 通用 N 维点类型（使用 Eigen 列向量） + * + * @tparam T 标量类型（double, float 等） + * @tparam Dim 维度 + * + * @ingroup foundation + */ template using Point = Eigen::Matrix; +/** + * @brief 通用 N 维向量类型（与 Point 相同的列向量表示） + * + * 虽然类型相同，但语义上区分几何点（位置）和向量（方向/位移）。 + * + * @tparam T 标量类型 + * @tparam Dim 维度 + * + * @ingroup foundation + */ template using Vector = Eigen::Matrix; +/// 3×3 矩阵模板 template using Matrix3 = Eigen::Matrix; +/// 4×4 矩阵模板 template using Matrix4 = Eigen::Matrix; +/// 双精度二维点 using Point2D = Point; +/// 双精度三维点 using Point3D = Point; +/// 单精度二维点 using Point2f = Point; +/// 单精度三维点 using Point3f = Point; + +/// 双精度齐次四维点（用于仿射变换） using Point4D = Eigen::Matrix; +/// 双精度二维向量 using Vector2D = Vector; +/// 双精度三维向量 using Vector3D = Vector; +/// 单精度二维向量 using Vector2f = Vector; +/// 单精度三维向量 using Vector3f = Vector; +/// 双精度 3×3 矩阵 using Matrix3D = Matrix3; +/// 双精度 4×4 矩阵 using Matrix4D = Matrix4; } // namespace vde::foundation + +/** + * @defgroup foundation 基础模块 + * + * 基础模块提供 VDE 引擎的核心基础设施： + * - 数学类型定义（点、向量、矩阵） + * - 精确几何谓词（自适应双精度 + GMP 可选后端） + * - 区间算术验证 + * - 分层容差管理 + * - 内存池分配 + * - 文件 I/O（STL、OBJ、PLY、glTF） + * - 格式化工具与序列化 + * + * @{ + */ + +/** @} */ // end of foundation group diff --git a/include/vde/foundation/memory_pool.h b/include/vde/foundation/memory_pool.h index 40f271a..72ce224 100644 --- a/include/vde/foundation/memory_pool.h +++ b/include/vde/foundation/memory_pool.h @@ -5,12 +5,52 @@ namespace vde::foundation { +/** + * @brief 固定块大小的内存池 + * + * 使用空闲链表实现 O(1) 的分配和释放。按固定大小的 chunk 预先分配内存， + * 适合频繁创建/销毁同类型小对象的场景（如半边网格的顶点和面）。 + * + * @tparam T 池中存储的对象类型 + * + * @note 内存池不调用对象的构造函数/析构函数，仅管理原始内存 + * @note 适合 POD 类型或通过 placement new 手动管理生命周期的类型 + * @warning 不在单个对象释放时归还内存给操作系统；仅在池析构时统一释放所有 chunk + * + * @code{.cpp} + * MemoryPool pool(4096); + * Vertex* v = pool.allocate(); // O(1) + * // ... 使用 v ... + * pool.deallocate(v); // O(1) + * @endcode + * + * @ingroup foundation + */ template class MemoryPool { public: + /** + * @brief 构造内存池 + * @param chunk_size 每个内存块可容纳的对象数（默认 1024） + */ explicit MemoryPool(size_t chunk_size = 1024) : chunk_size_(chunk_size) {} + + /** + * @brief 析构内存池，释放所有已分配的内存块 + */ ~MemoryPool() { for (auto* p : chunks_) ::operator delete(p); } + /** + * @brief 从池中分配 n 个 T 对象的内存 + * + * 当前实现忽略 n 参数，每次始终分配一个对象。 + * 若空闲链表为空，自动分配新的 chunk。 + * + * @param n 对象数量（当前未使用，保留用于 API 兼容） + * @return 指向未初始化内存的指针 + * + * @pre n 应小于 chunk_size_ + */ T* allocate(size_t n = 1) { if (free_list_ == nullptr) allocate_chunk(); T* p = static_cast(free_list_); @@ -18,12 +58,22 @@ public: return p; } + /** + * @brief 将内存归还到池中 + * + * 将对象内存插入空闲链表头部，不调用析构函数。 + * + * @param p 之前从 allocate() 获取的指针 + */ void deallocate(T* p, size_t = 1) { *reinterpret_cast(p) = free_list_; free_list_ = p; } private: + /** + * @brief 分配一个新的内存块并初始化空闲链表 + */ void allocate_chunk() { void* chunk = ::operator new(chunk_size_ * sizeof(T)); chunks_.push_back(chunk); diff --git a/include/vde/foundation/predicates.h b/include/vde/foundation/predicates.h index 4d56791..47789c8 100644 --- a/include/vde/foundation/predicates.h +++ b/include/vde/foundation/predicates.h @@ -6,21 +6,46 @@ namespace vde::foundation { -/// Unified predicate interface — selects GMP or adaptive double at compile time +/** + * @brief 统一谓词接口 — 编译时选择 GMP 或自适应双精度后端 + * + * 根据 VDE_USE_GMP 宏自动选择 GmpExact 或 AdaptiveDouble 作为后端。 + * 上层代码只需使用 Predicates，无需关心具体实现。 + * + * @ingroup foundation + */ struct Predicates { #ifdef VDE_USE_GMP + /// 使用 GMP 精确有理数算术（零误差，较慢） using Backend = exact::GmpExact; #else + /// 使用自适应双精度算术（Shewchuk 算法，默认） using Backend = exact::AdaptiveDouble; #endif + /** + * @brief 二维定向测试 + * @param a,b,c 三点 + * @return 定向结果 + */ static exact::Orientation orient_2d(const Point2D& a, const Point2D& b, const Point2D& c) { return Backend::orient_2d(a, b, c); } + /** + * @brief 三维定向测试 + * @param a,b,c,d 四点 + * @return 定向结果 + */ static exact::Orientation orient_3d(const Point3D& a, const Point3D& b, const Point3D& c, const Point3D& d) { return Backend::orient_3d(a, b, c, d); } + /** + * @brief 圆内测试 + * @param a,b,c 三角形三顶点 + * @param d 测试点 + * @return 圆测试结果 + */ static exact::CircleTest in_circle(const Point2D& a, const Point2D& b, const Point2D& c, const Point2D& d) { return Backend::in_circle(a, b, c, d); @@ -28,12 +53,50 @@ struct Predicates { }; // ── Convenience wrappers (preferred API) ── + +/** + * @brief 二维定向测试便捷封装 + * + * 使用 Predicates::orient_2d，无需指定后端。 + * 这是推荐的对外 API。 + * + * @param a,b,c 三个二维点 + * @return Orientation::Clockwise / CounterClockwise / Collinear + * + * @code{.cpp} + * using namespace vde::foundation; + * Point2D a(0,0), b(1,0), c(0,1); + * auto result = orient(a, b, c); // CounterClockwise + * @endcode + * + * @ingroup foundation + */ inline exact::Orientation orient(const Point2D& a, const Point2D& b, const Point2D& c) { return Predicates::orient_2d(a, b, c); } + +/** + * @brief 三维定向测试便捷封装 + * + * @param a,b,c 定义平面的三点 + * @param d 测试点 + * @return Orientation::Clockwise / CounterClockwise / Collinear + * + * @ingroup foundation + */ inline exact::Orientation orient(const Point3D& a, const Point3D& b, const Point3D& c, const Point3D& d) { return Predicates::orient_3d(a, b, c, d); } + +/** + * @brief 圆内测试便捷封装 + * + * @param a,b,c 三角形三顶点 + * @param d 测试点 + * @return CircleTest::Inside / On / Outside + * + * @ingroup foundation + */ inline exact::CircleTest in_circle(const Point2D& a, const Point2D& b, const Point2D& c, const Point2D& d) { return Predicates::in_circle(a, b, c, d); } diff --git a/include/vde/foundation/serializer.h b/include/vde/foundation/serializer.h index 167ac8b..e9e46f2 100644 --- a/include/vde/foundation/serializer.h +++ b/include/vde/foundation/serializer.h @@ -6,33 +6,113 @@ namespace vde::foundation { -// Versioned binary format for fast save/load -// Header: 8 bytes magic + 4 bytes version + 4 bytes flags - -constexpr uint64_t VDE_MAGIC = 0x56444547454F4D00ULL; // "VDEGEOM\0" -constexpr uint32_t VDE_FORMAT_VERSION = 1; - +/** + * @brief 二进制序列化格式的网格数据结构 + * + * 存储网格的扁平表示，顶点坐标连续排列，面索引连续排列。 + * 所有面均为三角形（face_valence 固定为 3）。 + * + * @ingroup foundation + */ struct SerializedMesh { + /// 顶点数量 uint32_t vertex_count; + /// 面数量 uint32_t face_count; - std::vector vertices; // xyz xyz ... (flat) - std::vector face_indices; // v0 v1 v2 ... - std::vector face_valence; // 3 for all triangles + /// 顶点坐标，扁平排列：x0 y0 z0 x1 y1 z1 ... + std::vector vertices; + /// 面顶点索引，扁平排列：v0_0 v0_1 v0_2 v1_0 v1_1 v1_2 ... + std::vector face_indices; + /// 每面的顶点数（当前固定为 3，即纯三角形网格） + std::vector face_valence; }; +/** + * @brief 二进制序列化器 + * + * 提供网格的快速二进制序列化/反序列化，用于 VDE 本地格式（VDE 二进制格式）。 + * + * 文件格式： + * - 8 字节魔数 "VDEGEOM\0" (0x56444547454F4D00) + * - 4 字节格式版本号 + * - 4 字节标志位 + * - 序列化的网格数据 + * + * @note 该格式专为 VDE 内部快速读写设计，不适用于跨平台数据交换 + * @note 使用 glTF/OBJ/STL 进行外部交换，此格式用于本地缓存 + * + * @ingroup foundation + */ class BinarySerializer { public: - /// Serialize mesh to binary buffer + /** + * @brief 序列化网格到二进制缓冲区 + * + * @param mesh 要序列化的半边网格 + * @return 压缩后的二进制数据（含文件头） + * + * @code{.cpp} + * auto data = BinarySerializer::serialize(mesh); + * // data 可直接写入文件或通过网络传输 + * @endcode + */ static std::vector serialize(const mesh::HalfedgeMesh& mesh); - /// Deserialize mesh from binary buffer + /** + * @brief 从二进制缓冲区反序列化网格 + * + * @param data 之前由 serialize() 生成的二进制数据 + * @return 恢复的半边网格 + * + * @throws std::runtime_error 若数据校验失败（魔数/版本不匹配） + * + * @code{.cpp} + * auto mesh = BinarySerializer::deserialize(data); + * @endcode + */ static mesh::HalfedgeMesh deserialize(const std::vector& data); - /// Write to file + /** + * @brief 序列化并写入文件 + * + * 等价于 serialize() + 文件写入，一步完成。 + * + * @param path 输出文件路径 + * @param mesh 要保存的半边网格 + * @return true 写入成功 + * @return false 文件创建失败 + * + * @see read_file + */ static bool write_file(const std::string& path, const mesh::HalfedgeMesh& mesh); - /// Read from file + /** + * @brief 从文件读取并反序列化 + * + * 等价于文件读取 + deserialize()，一步完成。 + * + * @param path VDE 二进制格式文件路径 + * @return 恢复的半边网格 + * + * @throws std::runtime_error 若文件无法打开或格式无效 + * + * @see write_file + */ static mesh::HalfedgeMesh read_file(const std::string& path); }; +/** + * @brief VDE 二进制文件魔数："VDEGEOM\0" + * + * 用于识别 VDE 二进制网格文件。 + */ +constexpr uint64_t VDE_MAGIC = 0x56444547454F4D00ULL; // "VDEGEOM\0" + +/** + * @brief VDE 二进制格式当前版本号 + * + * 用于向前兼容：读取时检查版本号，若不一致则拒绝或执行迁移。 + */ +constexpr uint32_t VDE_FORMAT_VERSION = 1; + } // namespace vde::foundation diff --git a/include/vde/foundation/tolerance.h b/include/vde/foundation/tolerance.h index c48a87c..317dd98 100644 --- a/include/vde/foundation/tolerance.h +++ b/include/vde/foundation/tolerance.h @@ -6,44 +6,136 @@ namespace vde::foundation { -/// 分层容差：支持全局 + 局部覆盖 +/** + * @brief 分层容差管理类 + * + * 提供四级容差（绝对、相对、角度、吸附），支持全局实例和局部作用域覆盖。 + * 用于数值比较、几何判定和点合并时的精度控制。 + * + * 四级容差含义： + * - **绝对容差 (absolute)**: 对接近零的值使用，直接比较差值的绝对值 + * - **相对容差 (relative)**: 对大值使用，比较差值 / max(|a|, |b|) + * - **角度容差 (angular)**: 角度/方向比较时的阈值（弧度） + * - **吸附容差 (snapping)**: 点合并/对齐时的容差，通常比判定容差宽松 + * + * 预设精度等级： + * - Modeling(): 1e-6 — 通用建模精度 + * - Fitting(): 1e-3 — 近似拟合 + * - Intersection(): 1e-8 — 精确求交 + * - Snapping(): 1e-4 — 点吸附 + * + * @ingroup foundation + */ class Tolerance { public: + /** + * @brief 构造容差对象 + * + * @param absolute 绝对容差，默认 1e-6 + * @param relative 相对容差，默认 1e-8 + * @param angular 角度容差（弧度），默认 1e-8 + * @param snapping 吸附容差，默认 1e-4 + */ Tolerance(double absolute = 1e-6, double relative = 1e-8, double angular = 1e-8, double snapping = 1e-4) : absolute_(absolute), relative_(relative), angular_(angular), snapping_(snapping) {} // ── 判定 ── + + /** + * @brief 判断两点是否相等（基于绝对容差） + * + * 使用欧氏距离与绝对容差比较。 + * + * @param a 第一个点（Eigen 矩阵） + * @param b 第二个点（Eigen 矩阵） + * @return true 若两点的欧氏距离 < absolute_ + */ bool points_equal(const Eigen::MatrixXd& a, const Eigen::MatrixXd& b) const { return (a - b).norm() < absolute_; } + /** + * @brief 判断值是否接近零 + * + * @param value 测试值 + * @return true 若 |value| < absolute_ + */ template bool is_zero(T value) const { return std::abs(value) < absolute_; } + /** + * @brief 混合容差相等判断 + * + * 综合使用绝对和相对容差： + * |a - b| < absolute + relative * max(|a|, |b|) + * + * @tparam T 数值类型 + * @param a 第一个值 + * @param b 第二个值 + * @return true 若两值在容差范围内相等 + * + * @code{.cpp} + * Tolerance tol(1e-6, 1e-8); + * tol.equals(1.0000001, 1.0); // true + * tol.equals(1000.0001, 1000.0); // true (相对项起作用) + * @endcode + */ template bool equals(T a, T b) const { return std::abs(a - b) < absolute_ + relative_ * std::max(std::abs(a), std::abs(b)); } + /// 获取绝对容差 double absolute() const { return absolute_; } + /// 获取相对容差 double relative() const { return relative_; } + /// 获取角度容差 double angular() const { return angular_; } + /// 获取吸附容差 double snapping() const { return snapping_; } // ── 分层容差：继承并收紧 ── + + /** + * @brief 收紧容差 + * + * 所有容差乘以 factor（< 1），用于在子操作中提高精度要求。 + * + * @param factor 收紧因子（默认 0.1） + * @return 收紧后的新容差对象 + * + * @code{.cpp} + * auto rough = Tolerance(Tolerance::Fitting()); + * auto precise = rough.tighten(0.01); // 拟合 → 建模级 + * @endcode + */ Tolerance tighten(double factor = 0.1) const { return Tolerance(absolute_ * factor, relative_ * factor, angular_ * factor, snapping_ * factor); } + /** + * @brief 放宽容差 + * + * 所有容差乘以 factor（> 1），用于在上层操作中降低精度要求。 + * + * @param factor 放宽因子（默认 10.0） + * @return 放宽后的新容差对象 + */ Tolerance relax(double factor = 10.0) const { return Tolerance(absolute_ * factor, relative_ * factor, angular_ * factor, snapping_ * factor); } - // ── 合并：取更严格的那个 ── + /** + * @brief 合并两个容差，取更严格的（更小的） + * + * @param a 第一个容差 + * @param b 第二个容差 + * @return 每项取 min 的新容差 + */ static Tolerance min(const Tolerance& a, const Tolerance& b) { return Tolerance(std::min(a.absolute_, b.absolute_), std::min(a.relative_, b.relative_), @@ -52,14 +144,35 @@ public: } // ── 全局实例 ── + + /** + * @brief 获取全局容差实例 + * + * 线程局部存储，默认值为建模精度 (1e-6)。 + * + * @return 全局容差的引用 + */ static Tolerance& global() { return global_; } + + /** + * @brief 设置全局容差 + * + * @param tol 新的全局容差 + * + * @note 推荐使用 ToleranceScope 进行临时修改，确保自动恢复 + */ static void set_global(const Tolerance& tol) { global_ = tol; } // ── 默认建模容差等级 ── - static constexpr double Modeling() { return 1e-6; } // 建模级 - static constexpr double Fitting() { return 1e-3; } // 拟合级 - static constexpr double Intersection() { return 1e-8; } // 求交级 - static constexpr double Snapping() { return 1e-4; } // 吸附级 + + /// 建模级精度 (1e-6)：通用几何建模 + static constexpr double Modeling() { return 1e-6; } + /// 拟合级精度 (1e-3)：曲面拟合、近似计算 + static constexpr double Fitting() { return 1e-3; } + /// 求交级精度 (1e-8)：精确求交、布尔运算 + static constexpr double Intersection() { return 1e-8; } + /// 吸附级精度 (1e-4)：点合并、顶点焊接 + static constexpr double Snapping() { return 1e-4; } private: double absolute_; @@ -69,13 +182,40 @@ private: static Tolerance global_; }; -/// 局部容差堆栈：支持作用域内的临时容差 +/** + * @brief 局部容差作用域（RAII） + * + * 构造时设置新的全局容差，析构时自动恢复为之前的值。 + * 用于在特定代码块内临时修改容差。 + * + * @code{.cpp} + * void precise_operation() { + * ToleranceScope scope(Tolerance::Intersection()); + * // 在此作用域内使用求交级精度 + * // ... 布尔运算 ... + * } // scope 析构，自动恢复 + * @endcode + * + * @ingroup foundation + */ class ToleranceScope { public: + /** + * @brief 进入新的容差作用域 + * + * 保存当前全局容差并设置新值。 + * + * @param tol 此作用域内使用的容差 + */ explicit ToleranceScope(const Tolerance& tol) : previous_(Tolerance::global()) { Tolerance::set_global(tol); } + + /** + * @brief 退出容差作用域，恢复之前的全局容差 + */ ~ToleranceScope() { Tolerance::set_global(previous_); } + private: Tolerance previous_; }; diff --git a/include/vde/mesh/alpha_shapes.h b/include/vde/mesh/alpha_shapes.h index b93b782..da23020 100644 --- a/include/vde/mesh/alpha_shapes.h +++ b/include/vde/mesh/alpha_shapes.h @@ -8,12 +8,45 @@ using core::Point2D; using core::Point3D; using core::Vector3D; -/// Alpha Shapes: extract surface from point cloud -/// For alpha → ∞, returns convex hull; for alpha → 0, returns all faces -/// @param alpha Radius threshold: faces with circumradius > alpha are removed +/** + * @brief Alpha Shapes 点云表面重建 + * + * 从 3D 点云的 Delaunay 四面体化中提取表面，通过参数 α 控制细节级别： + * - α → ∞：返回凸包（所有边界四面体面都保留） + * - α → 0：返回所有 Delaunay 面（非常详细，含内部结构） + * - 中间值：剔除外接球半径 > α 的四面体，保留细节适中的表面 + * + * 算法：对每个 Delaunay 四面体，检查其外接球半径。 + * 若半径 ≤ α，四面体保留；否则剔除。最终表面 = 保留四面体与 + * 被剔除四面体（或外部）之间的边界三角形面。 + * + * @param points 输入 3D 点云 + * @param alpha α 半径阈值，通常取平均点间距的 1~3 倍 + * @return TetrahedronMesh 保留的四面体网格 + * + * @code{.cpp} + * // 从点云重建表面 + * auto tets = alpha_shapes(point_cloud, 0.5); + * auto [verts, tris] = alpha_shapes_surface(tets); + * @endcode + * @see alpha_shapes_surface delaunay_3d + * @ingroup mesh + */ TetrahedronMesh alpha_shapes(const std::vector& points, double alpha); -/// Convert alpha shapes tet mesh to triangle surface +/** + * @brief 将 Alpha Shapes 四面体网格的表面提取为三角形网格 + * + * 遍历四面体网格，提取所有仅出现在一个四面体上的面（边界面）， + * 即为重建的表面三角形。 + * + * @param tet_mesh Alpha Shapes 输出的四面体网格 + * @return {顶点数组, 三角形索引数组} + * + * @note 表面三角形法向指向四面体外部（右手定则） + * @see alpha_shapes + * @ingroup mesh + */ std::pair, std::vector>> alpha_shapes_surface(const TetrahedronMesh& tet_mesh); diff --git a/include/vde/mesh/cdt_2d.h b/include/vde/mesh/cdt_2d.h index 3f89213..4ff0719 100644 --- a/include/vde/mesh/cdt_2d.h +++ b/include/vde/mesh/cdt_2d.h @@ -8,8 +8,29 @@ using core::Point2D; using core::Point3D; using core::Vector3D; -/// Constrained Delaunay Triangulation (CDT) in 2D -/// Ensures specified constraint edges appear in the triangulation +/** + * @brief 约束 Delaunay 三角化（CDT） + * + * 在标准 Delaunay 三角化基础上强制指定的约束边必须出现在三角化中。 + * 通过边翻转（edge flip）反复消除与约束边相交的三角形边，直到所有 + * 约束边都成为三角形边。 + * + * 约束边在输出三角化中一定存在（作为三角形边），但不保证满足 + * 空圆性质——靠近约束边的三角形可能违反 Delaunay 条件。 + * + * @param points 输入 2D 点集 + * @param constraints 约束边列表，每项为 {起点索引, 终点索引} + * @return DelaunayResult 含约束边的三角化结果 + * + * @note 约束边不能相交（否则行为未定义）；自相交约束需先分割交点 + * @code{.cpp} + * // 带孔洞的三角化：约束边标记孔洞边界 + * auto result = constrained_delaunay_2d(points, + * {{0,1}, {1,2}, {2,1}, {3,0}}); + * @endcode + * @see delaunay_2d + * @ingroup mesh + */ DelaunayResult constrained_delaunay_2d( const std::vector& points, const std::vector>& constraints); diff --git a/include/vde/mesh/delaunay_2d.h b/include/vde/mesh/delaunay_2d.h index db94412..e73e424 100644 --- a/include/vde/mesh/delaunay_2d.h +++ b/include/vde/mesh/delaunay_2d.h @@ -8,12 +8,40 @@ using core::Point2D; using core::Point3D; using core::Vector3D; +/** + * @brief 2D Delaunay 三角化结果 + * @ingroup mesh + */ struct DelaunayResult { - std::vector vertices; - std::vector> triangles; // indices into vertices + std::vector vertices; ///< 顶点坐标（可能含超级三角形的辅助点） + std::vector> triangles; ///< 三角形索引（CCW 顺序） }; -/// Bowyer-Watson incremental Delaunay triangulation, O(n²) worst, O(n log n) typical +/** + * @brief 2D Delaunay 三角化（Bowyer-Watson 增量算法） + * + * 构建点集的 Delaunay 三角化：任意三角形的外接圆内不含其他点。 + * 该性质使三角化的最小角最大化，避免狭长三角形。 + * + * 算法流程： + * 1. 创建包围所有输入点的超级三角形 + * 2. 逐点插入： + * a. 找出外接圆包含新点的所有三角形（"坏三角形"） + * b. 删除坏三角形，形成多边形空洞 + * c. 将新点与空洞边界各边连接形成新三角形 + * 3. 删除与超级三角形顶点相关的三角形 + * + * @param points 输入 2D 点集（至少 3 个不共线点） + * @return DelaunayResult 三角化结果 + * + * @note 时间复杂度：平均 O(n log n)，最坏 O(n²) + * @code{.cpp} + * auto result = delaunay_2d({{0,0}, {1,0}, {0.5, 0.5}, {0,1}, {1,1}}); + * // result.triangles 为 Delaunay 三角形索引 + * @endcode + * @see constrained_delaunay_2d + * @ingroup mesh + */ DelaunayResult delaunay_2d(const std::vector& points); } // namespace vde::mesh diff --git a/include/vde/mesh/delaunay_3d.h b/include/vde/mesh/delaunay_3d.h index 07eee0d..e08aabf 100644 --- a/include/vde/mesh/delaunay_3d.h +++ b/include/vde/mesh/delaunay_3d.h @@ -8,12 +8,35 @@ using core::Point2D; using core::Point3D; using core::Vector3D; +/** + * @brief 四面体网格结构 + * @ingroup mesh + */ struct TetrahedronMesh { - std::vector vertices; - std::vector> tetrahedra; // 4 vertex indices + std::vector vertices; ///< 顶点坐标 + std::vector> tetrahedra; ///< 四面体索引，每个 {i0,i1,i2,i3} }; -/// 3D Delaunay tetrahedralization (Bowyer-Watson) +/** + * @brief 3D Delaunay 四面体化（Bowyer-Watson 增量算法） + * + * 将 3D 点集剖分为 Delaunay 四面体网格：任意四面体的外接球内不含其他顶点。 + * 算法是 2D 版本的自然推广： + * 1. 构建超级四面体包含所有点 + * 2. 逐点插入：删除外接球包含新点的四面体 → 空洞 → 重连 + * 3. 去除超级四面体相关单元 + * + * @param points 输入 3D 点集（至少 4 个不共面点） + * @return TetrahedronMesh 四面体网格 + * + * @note 时间复杂度与 2D 类似：平均 O(n log n)，最坏 O(n²) + * @note 输出包括内部和边界四面体；通过 alpha_shapes 可提取表面 + * @code{.cpp} + * auto tet_mesh = delaunay_3d(point_cloud); + * @endcode + * @see alpha_shapes + * @ingroup mesh + */ TetrahedronMesh delaunay_3d(const std::vector& points); } // namespace vde::mesh diff --git a/include/vde/mesh/geodesic.h b/include/vde/mesh/geodesic.h index 3a77849..5266110 100644 --- a/include/vde/mesh/geodesic.h +++ b/include/vde/mesh/geodesic.h @@ -7,10 +7,32 @@ using core::Point2D; using core::Point3D; using core::Vector3D; -/// Geodesic distance from source vertices using the Heat Method (Crane et al. 2013) -/// Returns per-vertex distance from the nearest source -/// @param sources Source vertex indices -/// @param t Time step (default = average edge length squared) +/** + * @brief 热方法（Heat Method）计算测地距离 + * + * 基于 Crane et al. 2013 的算法，通过求解三个线性系统高效计算 + * 曲面上任意点到源点集的测地距离（Geodesic Distance）。 + * + * 算法三步： + * 1. 热扩散：求解 (I - t·L)·u = u₀，其中 u₀ 在源点为 1 其余为 0 + * 2. 梯度归一化：X = -∇u / |∇u| + * 3. 泊松重建：L·φ = div(X)，得到测地距离 φ + * + * 相比精确测地线算法（如 MMP），热方法更快（~线性时间）且易于实现， + * 但精度受时间步长 t 和网格分辨率影响。 + * + * @param mesh 输入三角网格 + * @param sources 源顶点索引列表（距离为 0 的顶点） + * @param t 时间步长，< 0 时自动设为平均边长平方 + * @return 每顶点到最近源点的测地距离 + * + * @note 要求网格为连通流形；时间步长 t 越小精度越高但数值越不稳定 + * @code{.cpp} + * auto dist = geodesic_distance(mesh, {0, 10, 20}); + * // dist[i] = 顶点 i 到 {0,10,20} 中最近顶点的测地距离 + * @endcode + * @ingroup mesh + */ std::vector geodesic_distance(const HalfedgeMesh& mesh, const std::vector& sources, double t = -1.0); diff --git a/include/vde/mesh/halfedge_mesh.h b/include/vde/mesh/halfedge_mesh.h index 87809f0..b542922 100644 --- a/include/vde/mesh/halfedge_mesh.h +++ b/include/vde/mesh/halfedge_mesh.h @@ -9,67 +9,243 @@ using core::AABB3D; using core::Vector3D; using core::Point3D; +/** + * @brief 半边数据结构元素 + * + * 每条半边存储：起点顶点索引、所属面索引、前驱/后继/对侧半边索引。 + * 对侧半边通过 opposite_index 配对，相邻两半边索引总是 {2k, 2k+1}。 + * + * @ingroup mesh + */ struct Halfedge { - int vertex_index = -1; - int face_index = -1; - int next_index = -1; - int prev_index = -1; - int opposite_index = -1; + int vertex_index = -1; ///< 半边起点（终点 = next 半边的起点） + int face_index = -1; ///< 所属面的索引（边界半边为 -1） + int next_index = -1; ///< 沿面的下一条半边 + int prev_index = -1; ///< 沿面的上一条半边 + int opposite_index = -1; ///< 对侧半边索引（反方向半边） }; +/** + * @brief 面元素 + * + * 存储面的起始半边索引和边数（阶）。 + * + * @ingroup mesh + */ struct Face { - int halfedge_index = -1; - int valence = 3; + int halfedge_index = -1; ///< 该面的任意一条半边 + int valence = 3; ///< 边数（三角面 = 3） }; +/** + * @brief 半边数据结构三角网格 + * + * 核心拓扑数据结构，支持 O(1) 访问面的邻面、顶点的邻面环、 + * 边界检测等拓扑查询。数据内部一致性由构造方法保证。 + * + * 数据成员： + * - vertices_: 顶点位置数组 + * - halfedges_: 半边数组，相邻索引为对侧半边（0↔1, 2↔3, ...） + * - faces_: 面数组，每个面对应一条起始半边 + * - face_normals_ / vertex_normals_: 延迟计算的法向缓存 + * + * @ingroup mesh + */ class HalfedgeMesh { public: HalfedgeMesh() = default; - // ── Construction ── + // ── Construction ────────────────────────────────────── + + /** + * @brief 清空所有数据 + */ void clear(); + + /** + * @brief 添加顶点 + * @param p 3D 坐标 + * @return 新顶点的索引（0-based） + */ int add_vertex(const Point3D& p); + + /** + * @brief 添加面 + * @param vertex_indices 按逆时针顺序排列的顶点索引（至少 3 个） + * @return 新面的索引，失败返回 -1 + * @note 自动创建/查找半边并设置对侧关系；面法向遵循右手定则 + */ int add_face(const std::vector& vertex_indices); + + /** + * @brief 从三角形列表批量构建半边形网格 + * @param verts 顶点位置数组 + * @param tris 三角形索引数组，每个 {i0,i1,i2} + * @note 等价于反复调用 add_vertex + add_face，但内部做批处理优化 + * @code{.cpp} + * HalfedgeMesh mesh; + * mesh.build_from_triangles(vertices, triangles); + * @endcode + */ void build_from_triangles(const std::vector& verts, const std::vector>& tris); - // ── Accessors ── + // ── Accessors ───────────────────────────────────────── + + /** + * @brief 顶点数量 + * @return vertices_ 大小 + */ [[nodiscard]] size_t num_vertices() const { return vertices_.size(); } + + /** + * @brief 面数量 + * @return faces_ 大小 + */ [[nodiscard]] size_t num_faces() const { return faces_.size(); } + + /** + * @brief 边数量 + * @return halfedges_.size() / 2 + */ [[nodiscard]] size_t num_edges() const { return halfedges_.size() / 2; } + + /** + * @brief 访问顶点 + * @param idx 顶点索引，0 ≤ idx < num_vertices() + * @return 顶点坐标只读引用 + */ [[nodiscard]] const Point3D& vertex(size_t idx) const { return vertices_[idx]; } + + /** + * @brief 访问面 + * @param idx 面索引 + * @return Face 只读引用 + */ [[nodiscard]] const Face& face(size_t idx) const { return faces_[idx]; } + + /** + * @brief 访问半边 + * @param idx 半边索引 + * @return Halfedge 只读引用 + */ [[nodiscard]] const Halfedge& halfedge(size_t idx) const { return halfedges_[idx]; } + + /** + * @brief 更新顶点位置 + * @param idx 顶点索引 + * @param p 新坐标 + * @note 标记法向缓存为脏，下次查询时重新计算 + */ void set_vertex(size_t idx, const Point3D& p) { vertices_[idx] = p; normals_dirty_ = true; } - // ── Topology queries ── + // ── Topology queries ────────────────────────────────── + + /** + * @brief 获取面的所有顶点索引（按环绕顺序） + * @param fi 面索引 + * @return 顶点索引序列 + */ [[nodiscard]] std::vector face_vertices(int fi) const; + + /** + * @brief 获取顶点的邻面索引环 + * @param vi 顶点索引 + * @return 所有以 vi 为顶点的面索引 + */ [[nodiscard]] std::vector vertex_faces(int vi) const; + + /** + * @brief 是否为边界半边 + * @param hei 半边索引 + * @return 没有所属面时为 true + */ [[nodiscard]] bool is_boundary_edge(int hei) const { return halfedges_[hei].face_index < 0; } + + /** + * @brief 是否为边界顶点 + * @param vi 顶点索引 + * @return 邻接任何边界半边时为 true + */ [[nodiscard]] bool is_boundary_vertex(int vi) const; - // ── Normals ── + // ── Normals ─────────────────────────────────────────── + + /** + * @brief 面法向（几何法向，非归一化？由实现决定） + * @param fi 面索引 + * @return (v1-v0)×(v2-v0) 归一化结果 + */ [[nodiscard]] Vector3D face_normal(int fi) const; + + /** + * @brief 顶点法向（邻面法向的面积加权平均） + * @param vi 顶点索引 + * @return 归一化顶点法向 + */ [[nodiscard]] Vector3D vertex_normal(int vi) const; + + /** + * @brief 重新计算所有面法向和顶点法向 + * @note 修改顶点位置后自动标记脏位，调用此方法触发重新计算 + */ void update_normals(); - // ── Bounds ── + // ── Bounds ──────────────────────────────────────────── + + /** + * @brief 网格包围盒 + * @return AABB3D 轴对齐包围盒 + */ [[nodiscard]] AABB3D bounds() const; - // ── Iterators for range-based for ── + // ── Iterators for range-based for ───────────────────── + + /** + * @brief 面迭代器范围（支持 for-range） + */ class FaceRange; + + /** + * @brief 顶点单邻环迭代器 + */ class VertexOneRing; + + /** + * @brief 获取所有面范围（range-based for 支持） + * @return FaceRange 对象 + * @code{.cpp} + * for (auto& f : mesh.faces_range()) { ... } + * @endcode + */ [[nodiscard]] FaceRange faces_range() const; + + /** + * @brief 获取顶点的 1-ring 邻域迭代器 + * @param vi 顶点索引 + * @return VertexOneRing 对象 + * @code{.cpp} + * for (int vj : mesh.vertex_one_ring(vi)) { + * // vj 是 vi 的直接邻居 + * } + * @endcode + */ [[nodiscard]] VertexOneRing vertex_one_ring(int vi) const; private: - std::vector vertices_; - std::vector halfedges_; - std::vector faces_; - std::vector face_normals_; - std::vector vertex_normals_; - bool normals_dirty_ = true; + std::vector vertices_; ///< 顶点位置 + std::vector halfedges_; ///< 半边数组（成对存储） + std::vector faces_; ///< 面数组 + std::vector face_normals_; ///< 面法向缓存 + std::vector vertex_normals_; ///< 顶点法向缓存 + bool normals_dirty_ = true; ///< 法向脏标记 + /** + * @brief 查找或创建半边 (v0→v1) + * @param v0 起点顶点索引 + * @param v1 终点顶点索引 + * @return 半边索引 + */ int find_or_create_edge(int v0, int v1); }; diff --git a/include/vde/mesh/marching_cubes.h b/include/vde/mesh/marching_cubes.h index 0c9da93..4dc6051 100644 --- a/include/vde/mesh/marching_cubes.h +++ b/include/vde/mesh/marching_cubes.h @@ -3,30 +3,82 @@ #include #include #include +#include +#include namespace vde::mesh { using core::Point2D; using core::Point3D; using core::Vector3D; +/** + * @brief Marching Cubes 输出网格 + * @ingroup mesh + */ struct MCMesh { - std::vector vertices; - std::vector> triangles; + std::vector vertices; ///< 顶点坐标 + std::vector> triangles; ///< 三角形索引 }; -/// Marching Cubes for scalar field f(x,y,z) -/// Extracts isosurface at value = iso_level within the bounding box +/** + * @brief Marching Cubes 等值面提取 + * + * 从标量场 f(x,y,z) 中提取等值面 S = { (x,y,z) | f(x,y,z) = iso_level }。 + * 将包围盒划分为 resolution×resolution×resolution 体素网格， + * 在每个体素的 8 个角求值 f，查表确定等值面穿越该体素的三角形配置。 + * + * @param f 标量场函数 f(x,y,z) → double + * @param iso_level 等值面值（默认 0） + * @param bmin 包围盒最小角 + * @param bmax 包围盒最大角 + * @param resolution 每轴采样分辨率（≥ 1），总网格 = resolution³ 个立方体 + * @return MCMesh 三角形网格 + * + * @note 使用经典的 256 查表法（Lorensen & Cline 1987） + * @code{.cpp} + * // 提取球面等值面 + * auto f = [](double x,double y,double z) { + * return sdf_sphere(x,y,z, 0,0,0, 1.0); + * }; + * auto mesh = marching_cubes(f, 0.0, + * {-2,-2,-2}, {2,2,2}, 64); + * @endcode + * @see sdf_sphere sdf_box + * @ingroup mesh + */ MCMesh marching_cubes(const std::function& f, double iso_level, const Point3D& bmin, const Point3D& bmax, int resolution); -/// SDF sphere +/** + * @brief SDF 球体函数 + * + * @param x,y,z 采样点坐标 + * @param cx,cy,cz 球心坐标 + * @param r 半径 + * @return 有符号距离：内部 < 0，表面 = 0，外部 > 0 + * + * @code{.cpp} + * auto sphere_sdf = [](double x,double y,double z) { + * return sdf_sphere(x,y,z, 0,0,0, 1.0); + * }; + * @endcode + * @ingroup mesh + */ inline double sdf_sphere(double x, double y, double z, double cx, double cy, double cz, double r) { return std::sqrt((x-cx)*(x-cx) + (y-cy)*(y-cy) + (z-cz)*(z-cz)) - r; } -/// SDF box +/** + * @brief SDF 立方体函数（轴对齐，中心在原点） + * + * @param x,y,z 采样点坐标 + * @param hx,hy,hz 半边长 + * @return 有符号距离 + * + * @ingroup mesh + */ inline double sdf_box(double x, double y, double z, double hx, double hy, double hz) { double dx = std::abs(x) - hx, dy = std::abs(y) - hy, dz = std::abs(z) - hz; return std::sqrt(std::max(dx,0.0)*std::max(dx,0.0) @@ -35,14 +87,34 @@ inline double sdf_box(double x, double y, double z, double hx, double hy, double + std::min(std::max({dx,dy,dz}), 0.0); } -/// SDF smooth union +/** + * @brief SDF 平滑并集（Smooth Union） + * + * 对两个 SDF 做平滑混合，k 控制过渡区宽度。 + * + * @param d1,d2 两个 SDF 值 + * @param k 平滑宽度（> 0），值越大过渡越平滑 + * @return 混合后的 SDF 值 + * + * @ingroup mesh + */ namespace { inline double lerp_impl(double a, double b, double t) { return a + t * (b - a); } } inline double sdf_smooth_union(double d1, double d2, double k) { double h = std::clamp(0.5 + 0.5*(d2-d1)/k, 0.0, 1.0); return lerp_impl(d2, d1, h) - k*h*(1.0-h); } -/// SDF smooth subtraction +/** + * @brief SDF 平滑差集（Smooth Subtraction） + * + * d1 减去 d2，过渡区由 k 控制。 + * + * @param d1,d2 两个 SDF 值 + * @param k 平滑宽度（> 0） + * @return 混合后的 SDF 值 + * + * @ingroup mesh + */ inline double sdf_smooth_subtraction(double d1, double d2, double k) { double h = std::clamp(0.5 - 0.5*(d2+d1)/k, 0.0, 1.0); return lerp_impl(d2, -d1, h) + k*h*(1.0-h); diff --git a/include/vde/mesh/mesh_boolean.h b/include/vde/mesh/mesh_boolean.h index 626bb03..5cec92c 100644 --- a/include/vde/mesh/mesh_boolean.h +++ b/include/vde/mesh/mesh_boolean.h @@ -6,8 +6,39 @@ using core::Point2D; using core::Point3D; using core::Vector3D; -enum class BooleanOp { Union, Intersection, Difference, SymDiff }; +/** + * @brief 布尔运算类型（CSG 组合） + * @ingroup mesh + */ +enum class BooleanOp { + Union, ///< A ∪ B 并集 + Intersection, ///< A ∩ B 交集 + Difference, ///< A \ B 差集 + SymDiff ///< A Δ B 对称差集 +}; +/** + * @brief 三角网格布尔运算 + * + * 对两个封闭三角网格执行 CSG（Constructive Solid Geometry）布尔操作。 + * 内部流程： + * 1. 计算两个网格的三角形-三角形交线 + * 2. 沿交线细分三角形（conforming 三角化） + * 3. 根据布尔类型对每个子面片做内/外分类 + * 4. 缝合属于结果的部分形成新网格 + * + * @param a 第一个三角网格（必须为封闭流形） + * @param b 第二个三角网格（必须为封闭流形） + * @param op 布尔运算类型 + * @return 结果三角网格 + * + * @note 要求输入网格为封闭流形（watertight），否则分类不可靠 + * @code{.cpp} + * auto result = mesh_boolean(cube, sphere, BooleanOp::Difference); + * // result = 立方体减去球体的交集部分 + * @endcode + * @ingroup mesh + */ HalfedgeMesh mesh_boolean(const HalfedgeMesh& a, const HalfedgeMesh& b, BooleanOp op); } // namespace vde::mesh diff --git a/include/vde/mesh/mesh_curvature.h b/include/vde/mesh/mesh_curvature.h index bac0b6d..e14bb1a 100644 --- a/include/vde/mesh/mesh_curvature.h +++ b/include/vde/mesh/mesh_curvature.h @@ -7,15 +7,42 @@ using core::Point2D; using core::Point3D; using core::Vector3D; +/** + * @brief 离散曲率计算结果 + * @ingroup mesh + */ struct CurvatureResult { - std::vector gaussian; // per-vertex Gaussian curvature - std::vector mean; // per-vertex mean curvature - std::vector area; // per-vertex Voronoi area + std::vector gaussian; ///< 每顶点高斯曲率 K = κ₁·κ₂ + std::vector mean; ///< 每顶点平均曲率 H = (κ₁+κ₂)/2 + std::vector area; ///< 每顶点 Voronoi 面积（曲率归一化权重） }; -/// Discrete curvature using Cotan formula (Meyer et al. 2003) -/// Gaussian: K(v) = (2π - Σθ_j) / A_v (interior vertices) -/// Mean: H(v) = ||Σ(cot α + cot β)·e|| / (2·A_v) +/** + * @brief 离散曲率计算（Cotan 公式，Meyer et al. 2003） + * + * 对三角网格每个顶点计算离散高斯曲率和平均曲率。 + * + * 高斯曲率（内点）： + * K(v) = (2π - Σ_j θ_j) / A_v + * 其中 θ_j 是顶点 v 处面角的总和，A_v 是 Voronoi 面积 + * + * 平均曲率： + * H(v) = || Σ_j (cot α_j + cot β_j) · e_j || / (2 · A_v) + * 其中 α_j,β_j 是与边 e_j 相对的两个角 + * + * 边界顶点使用角度缺损公式修正。 + * + * @param mesh 输入三角网格 + * @return CurvatureResult 每顶点曲率 + * + * @note 要求网格为流形；边界顶点曲率使用近似公式 + * @code{.cpp} + * auto crv = compute_curvature(mesh); + * double max_curvature = *std::max_element(crv.gaussian.begin(), crv.gaussian.end()); + * // 可用于特征检测、网格分割等 + * @endcode + * @ingroup mesh + */ CurvatureResult compute_curvature(const HalfedgeMesh& mesh); } // namespace vde::mesh diff --git a/include/vde/mesh/mesh_parameterize.h b/include/vde/mesh/mesh_parameterize.h index a8cd912..9cbeeea 100644 --- a/include/vde/mesh/mesh_parameterize.h +++ b/include/vde/mesh/mesh_parameterize.h @@ -8,12 +8,62 @@ using core::Point2D; using core::Point3D; using core::Vector3D; -/// Tutte embedding (harmonic parameterization) for a mesh with fixed boundary -/// Boundary vertices are mapped to a circle; interior vertices are solved via Laplacian -/// Returns per-vertex UV coordinates as Point2D +/** + * @brief Tutte 参数化（调和参数化） + * + * 将边界顶点固定映射到单位圆上，内部顶点通过求解 Laplace 方程 + * Δu = 0, Δv = 0 得到 UV 坐标。这是最经典、最稳定的网格参数化方法。 + * + * 算法： + * 1. 检测网格边界环 + * 2. 将边界顶点按弧长比例映射到单位圆上 + * 3. 对每个内部顶点求解均匀 Laplace（权重 w_ij = 1）： + * v_i = Σ_j v_j / deg(i) + * 4. 构造稀疏线性系统 LU 求解 + * + * 优点：极稳定，保证无翻转（对凸边界网格） + * 缺点：边界固定，内部可能扭曲（非共形） + * + * @param mesh 输入三角网格（需含边界） + * @return 每顶点 UV 坐标（Point2D），顺序与 mesh 顶点一致 + * + * @note 网格必须有边界（开网格）；闭曲面需先切割为拓扑圆盘 + * @code{.cpp} + * auto uv = tutte_parameterization(mesh); + * // uv[i] 为顶点 i 的 (u,v) 坐标，范围大致在 [-1,1]² + * @endcode + * @see lscm_parameterization + * @ingroup mesh + */ [[nodiscard]] std::vector tutte_parameterization(const HalfedgeMesh& mesh); -/// LSCM (Least Squares Conformal Maps) — simplified version +/** + * @brief LSCM（Least Squares Conformal Maps）参数化 + * + * 基于 Lévy et al. 2002 的最小二乘共形映射。通过最小化共形能量 + * 使参数化尽可能保角（角度畸变小），同时允许边界自由。 + * + * 与 Tutte 的区别： + * - 无需固定边界 → 边界自然展开，畸变更小 + * - 需要固定至少 2 个顶点消除刚体自由度（内部处理） + * - 共形性更好，适合纹理映射 + * - 不保证无翻转（局部极小可能翻转） + * + * 能量函数： + * E(u,v) = Σ_T || ∇(u+iv) ||² · A_T + * 对每个三角形 T 最小化梯度扰动平方和 + * + * @param mesh 输入三角网格 + * @return 每顶点 UV 坐标（Point2D） + * + * @note 简化实现，不包含完整 Lévy 论文中所有的退化处理 + * @code{.cpp} + * auto uv = lscm_parameterization(mesh); + * // 适合纹理 mapping + * @endcode + * @see tutte_parameterization + * @ingroup mesh + */ [[nodiscard]] std::vector lscm_parameterization(const HalfedgeMesh& mesh); } // namespace vde::mesh diff --git a/include/vde/mesh/mesh_quality.h b/include/vde/mesh/mesh_quality.h index 37fe0f0..8c2875f 100644 --- a/include/vde/mesh/mesh_quality.h +++ b/include/vde/mesh/mesh_quality.h @@ -3,93 +3,213 @@ #include #include #include +#include namespace vde::mesh { using core::Point2D; using core::Point3D; using core::Vector3D; -// ── Triangle quality ──────────────────────────────── +// ════════════════════════════════════════════════════════ +// Triangle quality +// ════════════════════════════════════════════════════════ +/** + * @brief 三角网格质量指标 + * @ingroup mesh + */ struct MeshQuality { - double min_angle_deg = 0; - double max_angle_deg = 0; - double avg_aspect_ratio = 0; - double max_aspect_ratio = 0; - size_t degenerate_faces = 0; + double min_angle_deg = 0; ///< 最小三角形内角（度），≥ 0 + double max_angle_deg = 0; ///< 最大三角形内角（度），≤ 180 + double avg_aspect_ratio = 0; ///< 平均长宽比（外接圆半径 / 最短边） + double max_aspect_ratio = 0; ///< 最大长宽比 + size_t degenerate_faces = 0; ///< 退化面数（面积 ≈ 0） }; +/** + * @brief 评估三角网格质量 + * + * 对网格中每个三角面计算内角、长宽比，统计全局最值及平均值。 + * 质量好的三角网格应有： + * - min_angle > 20°（CFD 模拟要求 > 25°） + * - max_angle < 120°（避免狭长退化） + * - avg_aspect_ratio 接近 1 + * + * @param mesh 输入三角网格 + * @return MeshQuality 质量指标汇总 + * + * @code{.cpp} + * auto q = evaluate_mesh_quality(mesh); + * if (q.min_angle_deg < 15.0) { + * // 需要 remesh 或优化 + * } + * @endcode + * @ingroup mesh + */ MeshQuality evaluate_mesh_quality(const HalfedgeMesh& mesh); -// ── Tetrahedron quality ───────────────────────────── +// ════════════════════════════════════════════════════════ +// Tetrahedron quality +// ════════════════════════════════════════════════════════ +/** + * @brief 四面体网格质量指标 + * @ingroup mesh + */ struct TetQuality { - double min_jacobian = 0.0; // scaled Jacobian (normalised, 0=b−ad … 1=regular) - double max_jacobian = 0.0; - double avg_jacobian = 0.0; - double min_dihedral_deg = 0.0; // minimum dihedral angle (degrees) - double max_dihedral_deg = 0.0; // maximum dihedral angle - double max_aspect_ratio = 0.0; // circumsphere-radius / shortest-edge - size_t degenerate_tets = 0; - size_t total_tets = 0; + double min_jacobian = 0.0; ///< 最小归一化 Jacobian [0, 1]，0=退化 1=正四面体 + double max_jacobian = 0.0; ///< 最大归一化 Jacobian + double avg_jacobian = 0.0; ///< 平均归一化 Jacobian + double min_dihedral_deg = 0.0; ///< 最小二面角（度），太小会导致刚度矩阵病态 + double max_dihedral_deg = 0.0; ///< 最大二面角（度），接近 180° 表示退化 + double max_aspect_ratio = 0.0; ///< 最大长宽比（外接球半径 / 最短边） + size_t degenerate_tets = 0; ///< 退化四面体数（体积 ≈ 0 或 Jacobian ≈ 0） + size_t total_tets = 0; ///< 四面体总数 }; -/// Evaluate tetrahedron quality from a triangle mesh. -/// This function interprets the mesh as a tetrahedral mesh where -/// each face belongs to a tet defined by the face vertices + face centroid -/// (suitable for closed triangle surfaces interpreted as volume boundaries). +/** + * @brief 评估三角网格作为四面体边界的质量 + * + * 将每个三角面与其面重心构造虚拟四面体，评估这些四面体的质量。 + * 适用于将封闭三角表面解释为体积边界时的质量检测。 + * + * @param mesh 封闭三角网格（volume boundary） + * @return TetQuality 四面体质量指标 + * + * @note 这不是真正的四面体网格质量评估——仅用于表面网格的代理指标 + * @ingroup mesh + */ TetQuality evaluate_tet_quality(const HalfedgeMesh& mesh); -// ── Generic element quality ───────────────────────── +// ════════════════════════════════════════════════════════ +// Generic element quality +// ════════════════════════════════════════════════════════ -enum class ElementType { Tri, Quad, Tet, Hex }; - -struct ElemQuality { - double min_jacobian = 0.0; - double max_jacobian = 0.0; - double avg_jacobian = 0.0; - size_t total_elements = 0; - size_t degenerate = 0; +/** + * @brief 单元类型枚举 + * @ingroup mesh + */ +enum class ElementType { + Tri, ///< 三角形 + Quad, ///< 四边形（Stub） + Tet, ///< 四面体 + Hex ///< 六面体（Stub） }; -/// Dispatch by element type; implemented for Tri and Tet (Quad/Hex return stub). +/** + * @brief 通用单元质量指标 + * @ingroup mesh + */ +struct ElemQuality { + double min_jacobian = 0.0; ///< 最小归一化 Jacobian + double max_jacobian = 0.0; ///< 最大归一化 Jacobian + double avg_jacobian = 0.0; ///< 平均归一化 Jacobian + size_t total_elements = 0; ///< 单元总数 + size_t degenerate = 0; ///< 退化单元数 +}; + +/** + * @brief 按单元类型分派的网格质量评估 + * + * @param mesh 输入网格 + * @param type 单元类型 + * @return ElemQuality 质量指标 + * + * @note 目前仅 Tri 和 Tet 完整实现；Quad / Hex 返回 stub（全零） + * @ingroup mesh + */ ElemQuality evaluate_element_quality(const HalfedgeMesh& mesh, ElementType type); -// ── Quality report ────────────────────────────────── +// ════════════════════════════════════════════════════════ +// Quality report +// ════════════════════════════════════════════════════════ +/** + * @brief 详细质量报告 + * @ingroup mesh + */ struct QualityReport { - std::string element_type; // "triangle" / "tetrahedron" - size_t total = 0; - size_t degenerate = 0; + std::string element_type; ///< "triangle" 或 "tetrahedron" - // Binned histogram (quality ∈ [0,1]) + size_t total = 0; ///< 单元总数 + size_t degenerate = 0; ///< 退化单元数 + + /// 质量直方图（k=10 个区间，[0,0.1), [0.1,0.2), ..., [0.9,1.0]） static constexpr int kBins = 10; std::array histogram{}; - // Grade counts (A=excellent B=good C=acceptable D=poor F=fail) + /// 分级计数（A=优秀 B=良好 C=合格 D=差 F=不合格） size_t grade_A = 0, grade_B = 0, grade_C = 0, grade_D = 0, grade_F = 0; - double min_quality = 0.0; - double max_quality = 0.0; - double avg_quality = 0.0; - double stddev_quality = 0.0; + double min_quality = 0.0; ///< 最低质量值 + double max_quality = 0.0; ///< 最高质量值 + double avg_quality = 0.0; ///< 平均质量值 + double stddev_quality = 0.0; ///< 质量标准差 }; -/// Generate a full quality report from a triangle surface mesh. +/** + * @brief 生成三角网格的完整质量报告 + * + * 包含直方图、分级统计及基本统计量。 + * 分级标准（按归一化 Jacobian）： + * - A: [0.9, 1.0] — 优秀 + * - B: [0.7, 0.9) — 良好 + * - C: [0.5, 0.7) — 合格 + * - D: [0.3, 0.5) — 差 + * - F: [0.0, 0.3) — 不合格 + * + * @param mesh 输入三角网格 + * @return QualityReport 完整质量报告 + * + * @ingroup mesh + */ QualityReport mesh_quality_report(const HalfedgeMesh& mesh); -/// Generate a full quality report from a tetrahedral mesh. -/// The mesh is treated as having 4-vertex faces (tets) in the -/// same index order as the surface-mesh face loop. +/** + * @brief 生成四面体网格的完整质量报告 + * + * 将表面网格的面解释为四面体单元的边界（面 + 面心 = 虚拟四面体）。 + * + * @param mesh 封闭三角网格 + * @return QualityReport 完整质量报告 + * + * @ingroup mesh + */ QualityReport mesh_quality_report_tet(const HalfedgeMesh& mesh); -// ── Utility ───────────────────────────────────────── +// ════════════════════════════════════════════════════════ +// Utility +// ════════════════════════════════════════════════════════ -/// Compute the scaled Jacobian (determinant / (product of edge-lengths)) -/// for a triangle; returns [−1,1], 1 = equilateral. +/** + * @brief 三角形归一化 Jacobian + * + * 计算三角形 (a,b,c) 的带符号面积，归一化到 [-1, 1]。 + * 值 1 = 等边三角形（最优），0 = 退化（共线），-1 = 翻转（法向反了）。 + * + * 公式：J = det([b-a, c-a, n]) / (|b-a|·|c-a|·|b-c|? or area-based?) + * 实际实现使用面积与理想等边三角形面积之比。 + * + * @param a,b,c 三角形三个顶点 + * @return 归一化 Jacobian [-1, 1] + * + * @ingroup mesh + */ double tri_scaled_jacobian(const Point3D& a, const Point3D& b, const Point3D& c); -/// Compute the scaled Jacobian for a tetrahedron; returns [0,1]. +/** + * @brief 四面体归一化 Jacobian + * + * 计算四面体 (a,b,c,d) 的体积 Jacobian，归一化到 [0, 1]。 + * 值 1 = 正四面体（最优），0 = 退化（共面）。 + * + * 公式：J = det([b-a, c-a, d-a]) / (归一化因子) + * + * @param a,b,c,d 四面体四个顶点 + * @return 归一化 Jacobian [0, 1] + * + * @ingroup mesh + */ double tet_scaled_jacobian(const Point3D& a, const Point3D& b, const Point3D& c, const Point3D& d); diff --git a/include/vde/mesh/mesh_repair.h b/include/vde/mesh/mesh_repair.h index 6097fb4..d0476ac 100644 --- a/include/vde/mesh/mesh_repair.h +++ b/include/vde/mesh/mesh_repair.h @@ -6,13 +6,40 @@ using core::Point2D; using core::Point3D; using core::Vector3D; +/** + * @brief 网格修复选项 + * @ingroup mesh + */ struct RepairOptions { - bool fill_holes = true; - bool remove_duplicates = true; - bool fix_orientation = true; - bool remove_degenerate = true; + bool fill_holes = true; ///< 是否填充孔洞（补面） + bool remove_duplicates = true; ///< 是否合并重合顶点（距离 < ε） + bool fix_orientation = true; ///< 是否统一面法向（外翻） + bool remove_degenerate = true; ///< 是否删除退化面（面积为 0 的三角形） }; +/** + * @brief 自动网格修复 + * + * 对输入网格执行一系列修复操作，输出水密的流形网格： + * - fill_holes: 检测边界环并用 fan 三角化填充 + * - remove_duplicates: 空间哈希合并间距 < 1e-8 的重复顶点 + * - fix_orientation: 从最大连通分量开始遍历传播一致的半边朝向 + * - remove_degenerate: 剔除面积 < ε 的退化三角形 + * + * @param mesh 输入网格（可能有缺陷） + * @param opts 修复选项（默认全部启用） + * @return 修复后的网格 + * + * @note 修复为启发式算法，不保证 100% 成功；孔洞过大或非流形边可能仍有残留问题 + * @code{.cpp} + * auto fixed = repair_mesh(broken_mesh); // 默认全修复 + * + * RepairOptions opts; + * opts.fill_holes = false; // 仅做去重 + 定向 + * auto cleaned = repair_mesh(mesh, opts); + * @endcode + * @ingroup mesh + */ HalfedgeMesh repair_mesh(const HalfedgeMesh& mesh, const RepairOptions& opts = {}); } // namespace vde::mesh diff --git a/include/vde/mesh/mesh_simplify.h b/include/vde/mesh/mesh_simplify.h index 95c8eea..29b6db0 100644 --- a/include/vde/mesh/mesh_simplify.h +++ b/include/vde/mesh/mesh_simplify.h @@ -6,14 +6,41 @@ using core::Point2D; using core::Point3D; using core::Vector3D; +/** + * @brief 网格简化选项 + * @ingroup mesh + */ struct SimplifyOptions { - double target_ratio = 0.5; // Target face ratio (0, 1) - int max_iterations = 100; - bool preserve_boundary = true; + double target_ratio = 0.5; ///< 目标面数比例 (0, 1]，0.5 表示保留 50% 的面 + int max_iterations = 100; ///< 最大迭代次数（边塌缩上限） + bool preserve_boundary = true; ///< 是否保护边界边不被塌缩 }; -/// QEM (Quadric Error Metrics) mesh simplification -/// Returns simplified mesh with ~target_ratio * original face count +/** + * @brief QEM（Quadric Error Metrics）网格简化 + * + * 基于 Garland & Heckbert 1997 的 QEM 算法，通过迭代边塌缩减少三角面数量。 + * 每条边维护一个 4×4 误差二次型 Q = Σ Q_f（面片的顶点-平面距离平方和）， + * 塌缩目标点为 argmin v^T Q v，塌缩代价为该最小值。 + * + * 算法流程： + * 1. 为每个顶点计算初始 Q（关联所有邻面的平面二次型之和） + * 2. 为每条边计算最优塌缩点及代价，插入最小堆 + * 3. 循环弹出最小代价边，执行塌缩，更新受影响边的代价 + * 4. 直到面数达到 target_ratio 或堆为空 + * + * @param mesh 输入三角网格 + * @param opts 简化选项 + * @return 简化后的网格 + * + * @note preserve_boundary 模式下，边界边通过惩罚项（大代价）被保护 + * @code{.cpp} + * SimplifyOptions opts; + * opts.target_ratio = 0.1; // 保留 10% 面 + * auto lod = simplify_mesh(high_res_mesh, opts); + * @endcode + * @ingroup mesh + */ HalfedgeMesh simplify_mesh(const HalfedgeMesh& mesh, const SimplifyOptions& opts = {}); } // namespace vde::mesh diff --git a/include/vde/mesh/mesh_smooth.h b/include/vde/mesh/mesh_smooth.h index b15fba0..7fc1ee7 100644 --- a/include/vde/mesh/mesh_smooth.h +++ b/include/vde/mesh/mesh_smooth.h @@ -6,38 +6,134 @@ using core::Point2D; using core::Point3D; using core::Vector3D; -enum class SmoothMethod { Laplacian, Taubin, HCLaplacian, Bilateral }; - -struct SmoothOptions { - int iterations = 10; - double lambda = 0.5; // Laplacian weight - double mu = -0.53; // Taubin shrink (mu < -lambda for Taubin) - - // HC Laplacian parameters - double hc_alpha = 0.0; // push-back strength (0 = default auto) - double hc_beta = 0.5; // push-forward blend - - // Bilateral parameters - double bilateral_sigma_c = 0.0; // spatial sigma (0 = auto from avg edge length) - double bilateral_sigma_n = 0.3; // normal sigma (radians) - int bilateral_iters = 5; - - SmoothMethod method = SmoothMethod::Laplacian; +/** + * @brief 网格光顺方法 + * @ingroup mesh + */ +enum class SmoothMethod { + Laplacian, ///< 标准拉普拉斯光顺：v ← v + λ·L(v)，快速但会收缩 + Taubin, ///< Taubin λ|μ 光顺：先正后负步交替，体积保持 + HCLaplacian, ///< HC 拉普拉斯（Humphrey's Classes）：推拉两步保持形状/体积 + Bilateral ///< 双边滤波：法向加权去噪，保持尖锐特征 }; -/// Generic dispatcher — delegates to specific smooth functions +/** + * @brief 网格光顺参数 + * @ingroup mesh + */ +struct SmoothOptions { + int iterations = 10; ///< 迭代次数 + double lambda = 0.5; ///< 拉普拉斯步长（Standard / Taubin 正步） + + /// Taubin 负步系数，典型值 mu = -(lambda + ε)，默认 -0.53 + /// 必须满足 mu < -lambda 以实现体积保持 + double mu = -0.53; + + /// HC 拉普拉斯参数 + double hc_alpha = 0.0; ///< 回推强度，0 = 自动计算（推荐） + double hc_beta = 0.5; ///< 前推混合系数 + + /// 双边滤波参数 + double bilateral_sigma_c = 0.0; ///< 空间 σ，0 = 自动按平均边长计算 + double bilateral_sigma_n = 0.3; ///< 法向 σ（弧度），控制特征保持的敏感度 + int bilateral_iters = 5; ///< 双边内部迭代次数 + + SmoothMethod method = SmoothMethod::Laplacian; ///< 光顺方法选择 +}; + +/** + * @brief 网格光顺（通用分发器） + * + * 根据 opts.method 分派到具体的光顺实现： + * - Laplacian → smooth_laplacian() + * - Taubin → smooth_taubin() + * - HCLaplacian → smooth_hc_laplacian() + * - Bilateral → smooth_bilateral() + * + * @param mesh 输入三角网格 + * @param opts 光顺选项 + * @return 光顺后的网格（顶点位置更新，拓扑不变） + * + * @code{.cpp} + * SmoothOptions opts; + * opts.method = SmoothMethod::Taubin; + * opts.iterations = 20; + * auto smooth = smooth_mesh(noisy_mesh, opts); + * @endcode + * @ingroup mesh + */ HalfedgeMesh smooth_mesh(const HalfedgeMesh& mesh, const SmoothOptions& opts = {}); -/// Standard Laplacian smoothing +/** + * @brief 标准拉普拉斯光顺 + * + * 均匀拉普拉斯：v ← v + λ·(avg(neighbor_positions) - v) + * 最简单但会引入体积收缩，适用于轻度平滑。 + * + * @param mesh 输入网格 + * @param iterations 迭代次数 + * @param lambda 步长系数 (0, 1] + * @return 光顺后的网格 + * + * @note 迭代过多会导致网格坍塌；对于体积关键的应用优先选 Taubin 或 HC + * @see smooth_taubin smooth_hc_laplacian + * @ingroup mesh + */ HalfedgeMesh smooth_laplacian(const HalfedgeMesh& mesh, int iterations, double lambda); -/// Taubin λ|μ smoothing (volume-preserving) +/** + * @brief Taubin λ|μ 光顺（体积保持） + * + * 两阶段交替： + * 1. 正步（收缩滤波）：v ← v + λ·L(v) + * 2. 负步（膨胀滤波）：v ← v + μ·L(v)，其中 μ < 0 且 |μ| > |λ| + * + * 传递函数在低频通带增益 ≈ 1，有效去噪同时保持体积。 + * + * @param mesh 输入网格 + * @param iterations 完整 λ/μ 循环次数 + * @param lambda 正步系数 (0, 1] + * @param mu 负步系数，需满足 mu < 0 且 |mu| > lambda + * @return 光顺后的网格 + * + * @note 典型值：λ = 0.5, μ = -0.53；|mu| 过大可能导致不稳定振荡 + * @ingroup mesh + */ HalfedgeMesh smooth_taubin(const HalfedgeMesh& mesh, int iterations, double lambda, double mu); -/// HC Laplacian (Humphrey's Classes) — two-pass with push-back to preserve volume +/** + * @brief HC 拉普拉斯光顺（Humphrey's Classes） + * + * 两遍处理以保持原始形状特征： + * 1. 前推：标准拉普拉斯光顺，记录位移 b_i = v_i' - v_i + * 2. 回推：v_i'' = v_i' - (α·b_i + β·avg(neighbor_b)) + * + * 相比标准拉普拉斯显著减轻体积收缩，同时比 Taubin 更好地保持尖锐特征。 + * + * @param mesh 输入网格 + * @param opts 光顺选项（使用 hc_alpha 和 hc_beta 字段） + * @return 光顺后的网格 + * + * @ingroup mesh + */ HalfedgeMesh smooth_hc_laplacian(const HalfedgeMesh& mesh, const SmoothOptions& opts = {}); -/// Bilateral mesh filtering — normal-weighted, preserves sharp features +/** + * @brief 双边网格滤波 + * + * 基于法向加权的各向异性去噪，保持尖锐边缘和特征： + * - 空间权重：w_c = exp(-d²/(2·σ_c²))，基于顶点间距 + * - 法向权重：w_n = exp(-θ²/(2·σ_n²))，基于法向夹角 + * - 更新：v ← v + Σ_j w_c(j)·w_n(j)·(v_j - v) / Σ_j w_c(j)·w_n(j) + * + * 相比于各向同性的拉普拉斯类方法，双边滤波在 CAD/机械零件网格上效果最优。 + * + * @param mesh 输入网格 + * @param opts 光顺选项（使用 bilateral_ 前缀字段） + * @return 光顺后的网格 + * + * @ingroup mesh + */ HalfedgeMesh smooth_bilateral(const HalfedgeMesh& mesh, const SmoothOptions& opts = {}); } // namespace vde::mesh diff --git a/include/vde/sdf/sdf_gradient.h b/include/vde/sdf/sdf_gradient.h index aeed226..a79c4d5 100644 --- a/include/vde/sdf/sdf_gradient.h +++ b/include/vde/sdf/sdf_gradient.h @@ -1,4 +1,28 @@ #pragma once +/** + * @file sdf_gradient.h + * @brief SDF 梯度计算（数值微分、解析梯度、链式法则） + * + * 提供两类梯度计算： + * + * 1. **空间梯度** ∇f(p) — SDF 在场中的方向导数（法向量方向） + * - 有限差分近似 + * - 各图元的解析梯度 + * - CSG 操作的链式法则 + * + * 2. **参数梯度** ∂f/∂θ — SDF 对形状参数的敏感度 + * - 有限差分计算 ∂f/∂radius, ∂f/∂height 等 + * - 用于基于梯度的形状优化 + * + * ## 梯度约定 + * + * 对于 SDF，梯度指向距离增加最快的方向。 + * 在表面上 (f=0)，∇f 等于外法向量（单位长度）。 + * 在外部 (f>0)，∇f 指向远离表面的方向。 + * 在内部 (f<0)，∇f 指向最近的表面。 + * + * @ingroup sdf + */ #include "vde/sdf/sdf_primitives.h" #include "vde/sdf/sdf_operations.h" #include "vde/core/point.h" @@ -14,11 +38,30 @@ using core::Vector3D; // Finite-Difference Gradient // ═══════════════════════════════════════════════ -/// Gradient of SDF at point p (normal direction) -/// Computed via central finite differences: (f(p+h) - f(p-h)) / (2h) -/// @param f SDF function (primitives or evaluate from tree) -/// @param p evaluation point -/// @param h step size (default 1e-6) +/** + * @brief SDF 的空间梯度（中心差分近似） + * + * 使用三轴中心有限差分: + * ∂f/∂x ≈ (f(x+h, y, z) - f(x-h, y, z)) / (2h) + * + * 二阶精度，适用于任意黑盒 SDF 函数。 + * + * @param f SDF 函数（接受 Point3D 参数） + * @param p 求值点 + * @param h 差分步长（默认 1e-6） + * @return 梯度向量 ∇f(p) + * + * @note 步长太小会导致舍入误差主导，太大导致截断误差主导。 + * 默认值 1e-6 在单精度场景下是较好的平衡。 + * @note 对于含不连续性的 SDF（如 sharp CSG），有限差分可能在边界产生错误梯度。 + * + * @code{.cpp} + * auto f = [](const Point3D& p) { return sphere(p, 1.0); }; + * Vector3D n = gradient(f, Point3D(1, 0, 0)); // ≈ (1, 0, 0) + * @endcode + * + * @see evaluate_with_gradient 一次性获取值和梯度 + */ [[nodiscard]] inline Vector3D gradient( const std::function& f, const Point3D& p, double h = 1e-6) @@ -35,13 +78,27 @@ using core::Vector3D; // Gradient Result Struct // ═══════════════════════════════════════════════ -/// Forward-mode AD: compute SDF value AND gradient in one pass -/// Returns {sdf_value, gradient_vector} +/** + * @brief 前向模梯度结果: SDF 值 + 梯度 + * + * 在一次计算中返回 SDF 值和空间梯度，避免重复求值。 + * 梯度向量在表面上即为单位外法向量的近似。 + */ struct GradResult { - double value; - Vector3D grad; // df/dx, df/dy, df/dz + double value; ///< SDF 值 f(p) + Vector3D grad; ///< 空间梯度 ∇f(p) = (df/dx, df/dy, df/dz) }; +/** + * @brief 同时计算 SDF 值和梯度的便捷函数 + * + * @param f SDF 函数 + * @param p 求值点 + * @param h 差分步长 + * @return value + grad 结构 + * + * @see gradient 仅计算梯度 + */ [[nodiscard]] inline GradResult evaluate_with_gradient( const std::function& f, const Point3D& p, double h = 1e-6) @@ -53,13 +110,36 @@ struct GradResult { // Parameter Gradients (Finite Difference) // ═══════════════════════════════════════════════ -/// Compute df/dradius for a sphere at point p +/** + * @brief 球体 SDF 对半径的参数梯度 ∂f/∂r + * + * 用中心差分计算: (sphere(p, r+h) - sphere(p, r-h)) / (2h) + * + * @param p 采样点 + * @param radius 球体半径 + * @param h 差分步长（默认 1e-6） + * @return ∂f/∂radius + * + * @note 对于球体，解析梯度为 -1（球面上），但此函数提供与实现无关的数值验证。 + * + * @see sphere_param_grad 封装版本 + */ [[nodiscard]] inline double sphere_radius_gradient(const Point3D& p, double radius, double h = 1e-6) { return (sphere(p, radius + h) - sphere(p, radius - h)) / (2.0 * h); } -/// Compute df/d(extents[i]) for a box (axis 0=x, 1=y, 2=z) +/** + * @brief 长方体 SDF 对各轴半尺寸的参数梯度 ∂f/∂extent[i] + * + * @param p 采样点 + * @param extents 半边长 + * @param axis 轴索引: 0=x, 1=y, 2=z + * @param h 差分步长 + * @return ∂f/∂extent[axis] + * + * @see box_param_grad 封装版本 + */ [[nodiscard]] inline double box_extent_gradient(const Point3D& p, const Point3D& extents, int axis, double h = 1e-6) { @@ -70,28 +150,52 @@ struct GradResult { return (box(p, ext_plus) - box(p, ext_minus)) / (2.0 * h); } -/// Compute df/d(major_radius) for a torus +/** + * @brief 圆环 SDF 对主半径的参数梯度 ∂f/∂major_r + * + * @param p 采样点 + * @param major_r 主半径 + * @param minor_r 副半径 + * @param h 差分步长 + * @return 参数梯度 + */ [[nodiscard]] inline double torus_major_gradient(const Point3D& p, double major_r, double minor_r, double h = 1e-6) { return (torus(p, major_r + h, minor_r) - torus(p, major_r - h, minor_r)) / (2.0 * h); } -/// Compute df/d(minor_radius) for a torus +/** + * @brief 圆环 SDF 对副半径的参数梯度 ∂f/∂minor_r + * + * @return 参数梯度 + */ [[nodiscard]] inline double torus_minor_gradient(const Point3D& p, double major_r, double minor_r, double h = 1e-6) { return (torus(p, major_r, minor_r + h) - torus(p, major_r, minor_r - h)) / (2.0 * h); } -/// Compute df/dheight for cylinder +/** + * @brief 圆柱 SDF 对高度的参数梯度 ∂f/∂height + * + * @param p 采样点 + * @param radius 圆柱半径 + * @param height 圆柱高度 + * @param h 差分步长 + * @return 参数梯度 + */ [[nodiscard]] inline double cylinder_height_gradient(const Point3D& p, double radius, double height, double h = 1e-6) { return (cylinder(p, radius, height + h) - cylinder(p, radius, height - h)) / (2.0 * h); } -/// Compute df/dradius for cylinder +/** + * @brief 圆柱 SDF 对半径的参数梯度 ∂f/∂radius + * + * @return 参数梯度 + */ [[nodiscard]] inline double cylinder_radius_gradient(const Point3D& p, double radius, double height, double h = 1e-6) { @@ -102,8 +206,20 @@ struct GradResult { // Analytical Gradients (Exact) // ═══════════════════════════════════════════════ -/// Exact gradient of sphere: radial unit vector (p - center) / |p - center| -/// For sphere centered at origin: p / |p| +/** + * @brief 球体 SDF 的解析梯度 + * + * **公式:** ∇f = p / |p| + * + * 即从球心指向采样点的单位径向向量。 + * 注意：球心处 p = (0,0,0)，梯度退化，返回零向量。 + * + * @param p 采样点 + * @param radius 球半径（当前未使用，保留以保持接口一致） + * @return 解析梯度（单位向量） + * + * @note 在表面上 (|p| = radius)，梯度即为外法向量。 + */ [[nodiscard]] inline Vector3D sphere_gradient_analytic(const Point3D& p, double /*radius*/) { double n = p.norm(); @@ -111,8 +227,19 @@ struct GradResult { return p / n; } -/// Exact gradient of box: sign(p) component-wise for exterior points -/// For interior points, the gradient is the unit normal of the nearest face +/** + * @brief 长方体 SDF 的解析梯度 + * + * **算法:** + * - 外部: 梯度分量非零当且仅当对应坐标超出 semi-extent，方向向外 + * - 内部: 梯度指向最近的表面（单分量） + * + * 注意：长方体 SDF 使用平滑近似，边缘处的梯度在内部/外部是分段常数。 + * + * @param p 采样点 + * @param half_extents 半边长 + * @return 解析梯度 + */ [[nodiscard]] inline Vector3D box_gradient_analytic(const Point3D& p, const Point3D& half_extents) { Vector3D g = Vector3D::Zero(); @@ -136,8 +263,20 @@ struct GradResult { return g; } -/// Exact gradient of torus -/// SDF: sqrt((sqrt(px^2 + pz^2) - major)^2 + py^2) - minor +/** + * @brief 圆环 SDF 的解析梯度 + * + * **SDF:** f = √((√(pₓ²+p𝓏²) - R_major)² + p_y²) - R_minor + * + * 梯度通过链式法则推导，将环向距离和 Y 方向分别求导。 + * + * **特殊情况:** 在 Y 轴上 (pₓ=p𝓏=0)，x 和 z 梯度分量退化为 0。 + * + * @param p 采样点 + * @param major_r 主半径 + * @param minor_r 副半径 + * @return 解析梯度 + */ [[nodiscard]] inline Vector3D torus_gradient_analytic(const Point3D& p, double major_r, double minor_r) { double rho = std::sqrt(p.x() * p.x() + p.z() * p.z()); @@ -158,7 +297,18 @@ struct GradResult { ); } -/// Exact gradient of cylinder (infinite along Y): normalize(p.x, 0, p.z) +/** + * @brief 无限圆柱（沿 Y 轴）的解析梯度 + * + * **公式:** ∇f = (pₓ, 0, p𝓏) / √(pₓ² + p𝓏²) + * + * 梯度仅含水平分量，方向为从轴线径向向外。 + * 在轴线上退化返回 (0,1,0)。 + * + * @param p 采样点 + * @param radius 圆柱半径（未使用，保留接口一致） + * @return 解析梯度 + */ [[nodiscard]] inline Vector3D cylinder_gradient_analytic(const Point3D& p, double /*radius*/) { double rho = std::sqrt(p.x() * p.x() + p.z() * p.z()); @@ -166,13 +316,41 @@ struct GradResult { return Vector3D(p.x() / rho, 0.0, p.z() / rho); } -/// Exact gradient of plane: normal +/** + * @brief 平面的解析梯度 + * + * **公式:** ∇f = normal + * + * 平面的梯度为常向量（法向量），与位置无关。 + * + * @param normal 单位法向量 + * @return 梯度（等于 normal） + * + * @note 这是最简单的解析梯度，因为平面 SDF 是线性的。 + */ [[nodiscard]] inline Vector3D plane_gradient_analytic(const Vector3D& normal) { return normal; } -/// Exact gradient of capsule: gradient of segment distance +/** + * @brief 胶囊体 SDF 的解析梯度 + * + * 梯度指向点 p 到最近线段上的投影方向。 + * + * **算法:** + * 1. 找到 p 在线段 ab 上的投影点 closest + * 2. 若 |p-closest| ≈ 0（点在轴上），返回零向量 + * 3. 否则梯度 = (p - closest) / |p - closest| + * + * 退化情况: 当 a ≈ b（线段退化为点）时，退化为球体的径向梯度。 + * + * @param p 采样点 + * @param a 线段起点 + * @param b 线段终点 + * @param radius 扫掠球半径 + * @return 解析梯度 + */ [[nodiscard]] inline Vector3D capsule_gradient_analytic(const Point3D& p, const Point3D& a, const Point3D& b, double radius) @@ -198,29 +376,79 @@ struct GradResult { // Chain Rule for CSG Operations // ═══════════════════════════════════════════════ -/// Gradient after union: gradient of the child that is closer +/** + * @brief 并集操作后的梯度传递 + * + * 取两个子对象中较近者的梯度（因为距离值来自较近者）。 + * + * @param grad_a 子对象 A 的梯度 + * @param grad_b 子对象 B 的梯度 + * @param d1 A 的 SDF 值 + * @param d2 B 的 SDF 值 + * @return 合并后的梯度 + * + * @warning 当 d1 ≈ d2 时，梯度可能在两个子对象间不连续跳变。 + * + * @see chain_smooth_union 避免不连续的平滑版本 + */ [[nodiscard]] inline Vector3D chain_union(const Vector3D& grad_a, const Vector3D& grad_b, double d1, double d2) { return (d1 < d2) ? grad_a : grad_b; } -/// Gradient after intersection: gradient of the child that is farther +/** + * @brief 交集操作后的梯度传递 + * + * 取两个子对象中较远者的梯度（因为距离值来自较远者）。 + * + * @return 合并后的梯度 + */ [[nodiscard]] inline Vector3D chain_intersection(const Vector3D& grad_a, const Vector3D& grad_b, double d1, double d2) { return (d1 > d2) ? grad_a : grad_b; } -/// Gradient after difference: max(d1, -d2) +/** + * @brief 差集操作后的梯度传递 + * + * 若 A 决定结果 (d1 > -d2)，使用 ∇A；否则使用 -∇B（B 的内部变成外部）。 + * + * @param grad_a 子对象 A 的梯度 + * @param grad_b 子对象 B 的梯度 + * @param d1 f_A + * @param d2 f_B + * @return 合并后的梯度 + */ [[nodiscard]] inline Vector3D chain_difference(const Vector3D& grad_a, const Vector3D& grad_b, double d1, double d2) { return (d1 > -d2) ? grad_a : Vector3D(-grad_b.x(), -grad_b.y(), -grad_b.z()); } -/// Gradient after smooth union with full chain rule -/// f = d1*(1-h) + d2*h - k*h*(1-h) where h = clamp(0.5+0.5*(d2-d1)/k, 0, 1) +/** + * @brief 平滑并集操作的梯度传递（含完整链式法则） + * + * 与普通并集的硬切换不同，本函数使用多项式混合函数的导数来平滑插值两个子对象的梯度。 + * + * **推导:** + * f = d₁(1-h) + d₂h - k·h(1-h)，h = clamp(0.5+0.5(d₂-d₁)/k, 0, 1) + * ∂f/∂d₁ = 2 - 3h + * ∂f/∂d₂ = 3h - 1 + * + * 在混合区域 (0 < h < 1) 内，权重平滑地从 A 过渡到 B。 + * + * @param grad_a 子对象 A 的梯度 + * @param grad_b 子对象 B 的梯度 + * @param d1 f_A + * @param d2 f_B + * @param k 混合参数 + * @return 平滑合并后的梯度 + * + * @see op_smooth_union 对应的平滑并集 SDF 操作 + * @see chain_union 硬切换版本 + */ [[nodiscard]] inline Vector3D chain_smooth_union(const Vector3D& grad_a, const Vector3D& grad_b, double d1, double d2, double k) { @@ -240,14 +468,33 @@ struct GradResult { // Chain Rule for Domain Transforms // ═══════════════════════════════════════════════ -/// Chain through translation: gradient unchanged (isometric) +/** + * @brief 平移变换后的梯度传递（梯度不变） + * + * 平移是等距变换（isometry），函数值不变， + * 因此子对象的梯度直接穿透。 + * + * @param child_result 变换空间中的子对象求值结果 + * @param offset 平移量（仅接口保留） + * @return 相同梯度 + */ [[nodiscard]] inline GradResult chain_translate(const GradResult& child_result, const Point3D& /*offset*/) { return child_result; } -/// Chain through rotation around Y axis (inverse rotation applied to gradient) -/// op_rotate applies R(-θ) to the point, so chain applies R(θ) to the gradient +/** + * @brief 绕 Y 轴旋转后的梯度传递 + * + * 由于 op_rotate 对采样点施加 R(-θ)，链式法则要求对梯度施加 R(θ)。 + * 即：∇f_world = R(θ) · ∇f_local + * + * @param child_result 旋转后坐标系中的子对象求值 + * @param angle_rad 旋转角度 + * @return 世界坐标系的梯度 + * + * @see op_rotate 对应的域变形操作 + */ [[nodiscard]] inline GradResult chain_rotate(const GradResult& child_result, double angle_rad) { GradResult result = child_result; @@ -262,7 +509,19 @@ struct GradResult { return result; } -/// Chain through non-uniform scale: divide gradient components by scale factors +/** + * @brief 非均匀缩放后的梯度传递 + * + * **公式:** ∇f_world = (∇f_child.x / sx, ∇f_child.y / sy, ∇f_child.z / sz) + * + * 因为 f_world(p) = f_child(p / s)，由链式法则得 ∂f_world/∂p = ∇f_child · diag(1/s)。 + * + * @note 缩放后梯度不再是单位长度，除非是均匀缩放。 + * + * @param child_result 缩放后坐标系中的子对象求值 + * @param factors 缩放因子 + * @return 世界坐标系的梯度 + */ [[nodiscard]] inline GradResult chain_scale(const GradResult& child_result, const Point3D& factors) { GradResult result = child_result; @@ -274,9 +533,28 @@ struct GradResult { return result; } -/// Chain through twist around Y axis -/// Twist: T(p) rotates (x,z) by +amount*p.y -/// Chain: J_T^T * child_grad where J_T includes y-dependence of the rotation angle +/** + * @brief 扭曲变换后的梯度传递（含完整雅可比） + * + * 扭曲 T(p) 绕 Y 轴旋转角度 amount·p_y。 + * 梯度传递需要对雅可比矩阵的转置做链式法则： + * + * **雅可比 J:** + * ``` + * ∂T_x/∂p_x = c, ∂T_x/∂p_y = -a·(p_x·s + p_z·c), ∂T_x/∂p_z = -s + * ∂T_y/∂p_x = 0, ∂T_y/∂p_y = 1, ∂T_y/∂p_z = 0 + * ∂T_z/∂p_x = s, ∂T_z/∂p_y = a·(p_x·c - p_z·s), ∂T_z/∂p_z = c + * ``` + * 其中 c = cos(a·p_y), s = sin(a·p_y), a = amount。 + * + * @param child_result 扭曲空间中的子对象求值 + * @param amount 单位高度的旋转量 + * @param p 原始采样点 + * @return 世界坐标系的梯度 + * + * @see op_twist 对应的域变形操作 + * @see chain_rotate 不含 Y 依赖的简化旋转梯度传递 + */ [[nodiscard]] inline GradResult chain_twist(const GradResult& child_result, double amount, const Point3D& p) { @@ -303,7 +581,19 @@ struct GradResult { return result; } -/// Chain through repeat: gradient unchanged (periodic domain) +/** + * @brief 周期重复后的梯度传递（梯度不变） + * + * 重复映射是平移变换的叠加，每个晶格内梯度相同（周期域）。 + * 在晶格边界处存在不连续性。 + * + * @param child_result 中心晶格内的子对象求值 + * @param cell 晶格尺寸（仅接口保留） + * @param p 原始采样点（仅接口保留） + * @return 相同梯度 + * + * @warning 在晶格边界处梯度不连续，可能影响基于梯度的优化。 + */ [[nodiscard]] inline GradResult chain_repeat(const GradResult& child_result, const Point3D& /*cell*/, const Point3D& /*p*/) { @@ -314,22 +604,41 @@ struct GradResult { // Parameter Gradient Vector (All Shape Params) // ═══════════════════════════════════════════════ -/// Parameter gradient struct: holds all possible shape-parameter sensitivities +/** + * @brief 参数梯度结构体 + * + * 存储 SDF 对所有形状参数的偏导数： + * ∂f/∂radius, ∂f/∂height, ∂f/∂extent[3], ∂f/∂major, ∂f/∂minor, + * ∂f/∂angle, ∂f/∂translate[3], ∂f/∂rotate, ∂f/∂scale[3] + * + * 用于基于梯度的形状优化和参数拟合。 + * + * @see sphere_param_grad, box_param_grad, cylinder_param_grad + */ struct ParamGrad { - double d_radius = 0.0; - double d_height = 0.0; - Point3D d_extents = Point3D::Zero(); - double d_major = 0.0; - double d_minor = 0.0; - double d_angle = 0.0; - Vector3D d_normal = Vector3D::Zero(); - double d_offset = 0.0; - Point3D d_translate = Point3D::Zero(); - double d_rotate = 0.0; - Point3D d_scale = Point3D::Zero(); + double d_radius = 0.0; ///< ∂f/∂radius + double d_height = 0.0; ///< ∂f/∂height + Point3D d_extents = Point3D::Zero(); ///< ∂f/∂(extent.x, extent.y, extent.z) + double d_major = 0.0; ///< ∂f/∂major_radius + double d_minor = 0.0; ///< ∂f/∂minor_radius + double d_angle = 0.0; ///< ∂f/∂angle_rad + Vector3D d_normal = Vector3D::Zero(); ///< ∂f/∂normal (无效分量清零) + double d_offset = 0.0; ///< ∂f/∂offset + Point3D d_translate = Point3D::Zero(); ///< ∂f/∂translate_offset + double d_rotate = 0.0; ///< ∂f/∂rotation_angle + Point3D d_scale = Point3D::Zero(); ///< ∂f/∂scale_factors }; -/// Compute parameter gradients for a sphere primitive +/** + * @brief 球体所有参数的梯度（当前仅 radius 有意义） + * + * @param p 采样点 + * @param radius 球体半径 + * @param h 差分步长 + * @return 参数梯度（仅 d_radius 非零） + * + * @see sphere_radius_gradient 底层数值差分 + */ [[nodiscard]] inline ParamGrad sphere_param_grad(const Point3D& p, double radius, double h = 1e-6) { ParamGrad pg; @@ -337,7 +646,18 @@ struct ParamGrad { return pg; } -/// Compute parameter gradients for a box primitive +/** + * @brief 长方体所有参数的梯度 + * + * 分别计算三个轴的半尺寸对 SDF 值的影响。 + * + * @param p 采样点 + * @param half_extents 半边长 + * @param h 差分步长 + * @return 参数梯度（d_extents 分量填充） + * + * @see box_extent_gradient 底层数值差分 + */ [[nodiscard]] inline ParamGrad box_param_grad(const Point3D& p, const Point3D& half_extents, double h = 1e-6) { @@ -348,7 +668,17 @@ struct ParamGrad { return pg; } -/// Compute parameter gradients for a cylinder primitive +/** + * @brief 圆柱所有参数的梯度 + * + * 分别计算半径和高度对 SDF 值的影响。 + * + * @param p 采样点 + * @param radius 圆柱半径 + * @param height 圆柱高度 + * @param h 差分步长 + * @return 参数梯度（d_radius, d_height 填充） + */ [[nodiscard]] inline ParamGrad cylinder_param_grad(const Point3D& p, double radius, double height, double h = 1e-6) { diff --git a/include/vde/sdf/sdf_operations.h b/include/vde/sdf/sdf_operations.h index 7f5d8e2..b06c892 100644 --- a/include/vde/sdf/sdf_operations.h +++ b/include/vde/sdf/sdf_operations.h @@ -1,4 +1,19 @@ #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 #include @@ -10,30 +25,150 @@ 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); @@ -42,17 +177,77 @@ using core::Vector3D; // ── 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) ── -/// Infinite repetition in a cell +/** + * @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; @@ -62,53 +257,232 @@ using core::Vector3D; wrap(p.z(), cell.z())); } -/// Mirror across X=0 plane, with optional offset +/** + * @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()); } -/// Mirror across Y=0 +/** + * @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()); } -/// Mirror across Z=0 +/** + * @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())); } -/// Translate space (inverse for SDF: move the world, not the object) +/** + * @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; } -/// Rotate around Y axis by angle_rad (inverse for SDF) +/** + * @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); } -/// Non-uniform scale (inverse for SDF) +/** + * @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()); } -/// Twist space around Y axis +/** + * @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); } -/// Bend space +/** + * @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()); } -/// Elongate (stretch space, eliminating interior along major axes) +/** + * @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( @@ -118,14 +492,56 @@ using core::Vector3D; ); } -/// Cheap bend (simplified bending) +/** + * @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); } -/// Displace distance field with sinusoidal noise +/** + * @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) * diff --git a/include/vde/sdf/sdf_primitives.h b/include/vde/sdf/sdf_primitives.h index b852076..1bd1b02 100644 --- a/include/vde/sdf/sdf_primitives.h +++ b/include/vde/sdf/sdf_primitives.h @@ -1,4 +1,29 @@ #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 #include @@ -11,21 +36,93 @@ using core::Vector3D; // Basic Primitives // ═══════════════════════════════════════════════ -/// Signed distance to sphere centered at origin +/** + * @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; } -/// Signed distance to axis-aligned box centered at origin -/// half_extents = (hx, hy, hz) — half-dimensions in each axis +/** + * @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); } -/// Signed distance to rounded box centered at origin -/// Chamfer radius r is subtracted from the box distance +/** + * @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() @@ -33,16 +130,62 @@ using core::Vector3D; - r; } -/// Signed distance to torus in XZ plane, Y-axis symmetry -/// major_radius = distance from origin to tube center -/// minor_radius = tube thickness radius +/** + * @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; } -/// Signed distance to capsule (line segment swept with sphere of given radius) -/// a, b = segment endpoints +/** + * @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; @@ -55,8 +198,32 @@ using core::Vector3D; return (pa - ba * h).norm() - radius; } -/// Signed distance to capped cylinder along Y axis, centered at origin -/// radius = cylinder radius, height = full height (from -h/2 to +h/2) +/** + * @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; @@ -65,15 +232,65 @@ using core::Vector3D; + std::max(d_y, 0.0) * std::max(d_y, 0.0)); } -/// Signed distance to plane through origin with given normal -/// Normal must be unit length; offset shifts plane along normal +/** + * @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; } -/// Signed distance to ellipsoid centered at origin -/// radii = semi-axis lengths (rx, ry, rz) -/// Uses the standard bounded-gradient approximation: (|p/radii| - 1) * min(radii) +/** + * @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(); @@ -91,9 +308,28 @@ using core::Vector3D; // Parametric Primitives // ═══════════════════════════════════════════════ -/// Signed distance to regular hexagonal prism in XZ plane, extruded along Y -/// radius = circumradius of hexagon (center to vertex) -/// height = full extrusion height (±h/2 along Y) +/** + * @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()); @@ -103,15 +339,57 @@ using core::Vector3D; return std::max(d_hex, std::abs(p.y()) - height * 0.5); } -/// Signed distance to infinite cylinder through origin along given axis -/// Axis must be unit length +/** + * @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; } -/// Signed distance to wedge — half of a box cut diagonally in XZ plane -/// w = width (X), h = height (Y), d = depth (Z) -/// The wedge occupies the half where z ≥ x inside the box +/** + * @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() @@ -125,21 +403,79 @@ using core::Vector3D; // 2D Extrusions & Revolution // ═══════════════════════════════════════════════ -/// Extrude a 2D SDF infinitely along Z axis -/// sdf_2d = pre-computed 2D SDF value at (p.x, p.y) +/** + * @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; } -/// Bounded extrusion: 2D SDF intersected with Z slab ±half_height -/// sdf_2d = pre-computed 2D SDF value at (p.x, p.y) +/** + * @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); } -/// Revolve a 2D profile around Y axis -/// sdf_2d = pre-computed 2D SDF value at (|p.xz|-offset, p.y) -/// offset = radial offset from Y axis +/** + * @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; } @@ -148,24 +484,94 @@ using core::Vector3D; // Complex Primitives (implemented in src/sdf/) // ═══════════════════════════════════════════════ -/// Signed distance to capped cone, centered at origin along Y axis -/// Apex at y = +height/2, base at y = -height/2 -/// angle_rad = half-angle at apex, height = full height +/** + * @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); -/// Signed distance to triangular prism with arbitrary triangle cross-section -/// Triangle defined by vertices a,b,c in XY plane, extruded ±height/2 along Z +/** + * @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); -/// Signed distance to link (two parallel toruses connected along X) -/// length = spacing between torus centers along X -/// major_r = major radius of each torus -/// minor_r = minor (tube) radius of each torus +/** + * @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); -/// Signed distance to infinite cone from apex along axis -/// axis must be unit length, angle_rad = half-angle +/** + * @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); diff --git a/include/vde/sdf/sdf_to_mesh.h b/include/vde/sdf/sdf_to_mesh.h index 2de6b51..73951c6 100644 --- a/include/vde/sdf/sdf_to_mesh.h +++ b/include/vde/sdf/sdf_to_mesh.h @@ -1,26 +1,91 @@ #pragma once +/** + * @file sdf_to_mesh.h + * @brief SDF 到三角网格的转换 + * + * 使用 Marching Cubes 算法将隐式曲面（SDF）转换为显式三角网格。 + * + * ## Marching Cubes 算法 + * + * 在指定的包围盒内以固定分辨率采样 SDF 场，通过查表确定每个体素内的等值面三角形配置。 + * + * **流程:** + * 1. 在 `[bmin, bmax]` 范围内构建分辨率 `res³` 的三维规则网格 + * 2. 在每个网格点求值 SDF + * 3. 利用 Marching Cubes 查找表在等值面处生成三角形 + * 4. 返回包含顶点和三角形索引的 MCMesh + * + * ## 注意事项 + * + * - 分辨率越高网格越精细，但计算量按 O(res³) 增长 + * - 推荐分辨率 64~256，取决于几何复杂度 + * - 结果网格需在 Marching Cubes 模块中定义（vde::mesh::MCMesh） + * + * @ingroup sdf + */ #include "vde/sdf/sdf_tree.h" #include "vde/mesh/marching_cubes.h" #include namespace vde::sdf { -/// Convert an SDF expression tree to a triangle mesh using marching cubes -/// @param root SDF expression tree root node -/// @param resolution Grid resolution per axis (e.g. 64 → 64³ grid) -/// @param iso_level Isosurface level (0 = surface, default) -/// @return Mesh with vertices and triangle indices +/** + * @brief 将 SDF 表达式树转换为三角网格 + * + * 使用 Marching Cubes 算法在 SDF 树的包围盒内采样并生成网格。 + * 包围盒由 estimate_bbox() 自动确定。 + * + * **算法步骤:** + * 1. 调用 estimate_bbox(root) 确定采样范围 + * 2. 构建分辨率 `resolution³` 体素网格 + * 3. 在每个体素角点求值 SDF + * 4. 用 Marching Cubes 提取 iso_level 等值面 + * + * @param root SDF 表达式树的根节点 + * @param resolution 每轴的网格分辨率（例如 64 表示 64³ 体素网格） + * @param iso_level 等值面水平（0.0 = 隐式表面） + * @return 包含顶点和三角形索引的 MCMesh + * + * @note 分辨率推荐为 2 的幂（32, 64, 128, 256）以获得最佳缓存性能。 + * @note 对于大尺寸但细节少的几何体，使用较低的 resolution 即可获得光滑结果。 + * + * @code{.cpp} + * auto root = SdfNode::sphere(2.0); + * auto mesh = sdf_to_mesh(root, 128); // 128³ 体素网格 + * // mesh.vertices 和 mesh.indices 可直接用于渲染或导出 + * @endcode + * + * @see sdf_to_mesh_lambda Lambda 版本的 SDF 转网格（Python 友好） + * @see estimate_bbox 包围盒估算 + */ [[nodiscard]] mesh::MCMesh sdf_to_mesh(const SdfNodePtr& root, int resolution = 64, double iso_level = 0.0); -/// Convert a lambda SDF f(x,y,z) to a triangle mesh -/// Convenience overload for direct function binding (useful from Python) -/// @param f SDF function f(x,y,z) → signed distance -/// @param bmin Lower corner of bounding box -/// @param bmax Upper corner of bounding box -/// @param resolution Grid resolution per axis -/// @param iso_level Isosurface level +/** + * @brief 将 Lambda SDF 函数转换为三角网格 + * + * 接受 C++ lambda 或 std::function 作为 SDF 定义，适合 Python 绑定场景。 + * 与 sdf_to_mesh() 不同，本函数需要调用方显式提供包围盒。 + * + * @param f SDF 函数，签名为 double(double x, double y, double z) → 有符号距离 + * @param bmin 采样包围盒最小角点 + * @param bmax 采样包围盒最大角点 + * @param resolution 每轴分辨率 + * @param iso_level 等值面水平（默认 0.0） + * @return MCMesh 三角网格 + * + * @code{.cpp} + * // 球体 lambda + * auto sphere_lambda = [](double x, double y, double z) { + * return std::sqrt(x*x + y*y + z*z) - 1.0; + * }; + * auto mesh = sdf_to_mesh_lambda(sphere_lambda, + * Point3D(-1.5, -1.5, -1.5), Point3D(1.5, 1.5, 1.5), 64); + * @endcode + * + * @see sdf_to_mesh 基于 SDF 树的网格转换 + */ [[nodiscard]] mesh::MCMesh sdf_to_mesh_lambda( const std::function& f, const Point3D& bmin, const Point3D& bmax, diff --git a/include/vde/sdf/sdf_tree.h b/include/vde/sdf/sdf_tree.h index b2f9b3a..4f0ed6b 100644 --- a/include/vde/sdf/sdf_tree.h +++ b/include/vde/sdf/sdf_tree.h @@ -1,4 +1,37 @@ #pragma once +/** + * @file sdf_tree.h + * @brief SDF 表达式树结构 + * + * 将 SDF 原语和 CSG 操作以树形数据结构表示，支持树遍历、求值和包围盒估算。 + * + * ## 架构层级 + * + * ``` + * SdfNode (树节点) — 单节点，含操作类型+参数+子节点 + * ├─ enum SdfOp — 操作枚举（图元/CSG/修饰/变形） + * ├─ struct SdfParams — 参数联合体 + * ├─ 静态工厂方法 — 创建各类型节点 + * └─ visit() — 前序遍历 + * + * evaluate() — 递归求值 + * estimate_bounds() / estimate_bbox() — 包围盒估算 + * ``` + * + * ## 使用模式 + * + * 通过静态工厂方法构建节点树，然后调用 evaluate() 进行求值： + * + * @code{.cpp} + * auto root = SdfNode::op_union( + * SdfNode::sphere(1.0), + * SdfNode::box(Point3D(0.5, 0.5, 0.5)) + * ); + * double d = evaluate(root, Point3D(1.0, 0.0, 0.0)); + * @endcode + * + * @ingroup sdf + */ #include "vde/sdf/sdf_primitives.h" #include "vde/sdf/sdf_operations.h" #include "vde/core/point.h" @@ -14,7 +47,16 @@ using core::Vector3D; class SdfNode; using SdfNodePtr = std::shared_ptr; -/// All node types in the CSG tree +/** + * @brief CSG 树的所有节点类型 + * + * 枚举涵盖三类操作： + * - **图元**（Sphere, Box, ...） — 叶子节点，产生基本几何体 + * - **CSG 布尔**（Union, Intersection, Difference） — 二叉树组合 + * - **修饰/变形**（Round, Translate, Twist, ...） — 单子节点变换 + * + * @see SdfParams 对应的参数结构 + */ enum class SdfOp : uint8_t { // Primitives Sphere, Box, RoundBox, Torus, Capsule, Cylinder, Cone, Plane, @@ -31,87 +73,208 @@ enum class SdfOp : uint8_t { Displace }; -/// Parameter payload — one struct to rule them all +/** + * @brief SDF 节点参数联合体 + * + * 一个结构体囊括所有可能的参数。不同操作类型使用不同字段子集： + * + * | 操作类型 | 使用字段 | + * |-------------------|--------------------------------------------| + * | Sphere | radius | + * | Box / RoundBox | extents | + * | Cylinder / Cone | radius, height, angle_rad | + * | Plane | normal, offset | + * | Torus | major_radius, minor_radius | + * | Capsule | pt_a, pt_b, radius | + * | Smooth CSG | blend_k | + * | Repeat | repeat_cell | + * | Translate | translate_offset | + * | Scale | scale_factors | + * | Displace | amplitude, frequency | + * | Twist/Bend | amount | + * | Onion | thickness | + * | Round | offset | + * + * @warning 不使用的字段保持默认值即可，但不要依赖未使用字段的默认值来做逻辑判断。 + */ struct SdfParams { - // Primitives - double radius = 1.0; - Point3D extents = Point3D(1, 1, 1); - Point3D pt_a = Point3D(0, -1, 0); - Point3D pt_b = Point3D(0, 1, 0); - double major_radius = 1.0; - double minor_radius = 0.3; - double angle_rad = 0.5; - double height = 2.0; - double thickness = 0.1; - Vector3D normal = Vector3D(0, 1, 0); - double offset = 0.0; - // Operations - double blend_k = 0.2; - double amount = 0.5; - double amplitude = 0.1; - double frequency = 1.0; - Point3D repeat_cell = Point3D(2, 2, 2); - Point3D translate_offset = Point3D(0, 0, 0); - Point3D scale_factors = Point3D(1, 1, 1); + // ── Primitive parameters ── + double radius = 1.0; ///< 球体/圆柱/胶囊/圆角半径 + Point3D extents = Point3D(1, 1, 1); ///< 盒子的半边长 + Point3D pt_a = Point3D(0, -1, 0); ///< 胶囊端点 A + Point3D pt_b = Point3D(0, 1, 0); ///< 胶囊端点 B + double major_radius = 1.0; ///< 圆环主半径 + double minor_radius = 0.3; ///< 圆环副半径 + double angle_rad = 0.5; ///< 锥体半角（弧度） + double height = 2.0; ///< 柱体/锥体/挤出高度 + double thickness = 0.1; ///< 壳层厚度（Onion） + Vector3D normal = Vector3D(0, 1, 0); ///< 平面法向量 + double offset = 0.0; ///< 平面偏移 / 通用偏移 + + // ── Operation parameters ── + double blend_k = 0.2; ///< 平滑 CSG 混合参数 + double amount = 0.5; ///< 扭曲/弯曲量 + double amplitude = 0.1; ///< 位移扰动幅值 + double frequency = 1.0; ///< 位移扰动频率 + Point3D repeat_cell = Point3D(2, 2, 2); ///< 重复晶格尺寸 + Point3D translate_offset = Point3D(0, 0, 0); ///< 平移量 + Point3D scale_factors = Point3D(1, 1, 1); ///< 缩放因子 }; -/// A single node in the SDF expression tree +/** + * @brief SDF 表达式树节点 + * + * 每个节点包含： + * - `op`: 操作类型 + * - `params`: 参数值 + * - `children`: 子节点列表（CSG 布尔 = 2 个子节点，修饰/变形 = 1 个，图元 = 0 个） + * - `name`: 可选名称（调试用） + * + * 使用 make_shared(op) 或静态工厂方法创建节点。 + * 所有工厂方法返回 SdfNodePtr（shared_ptr），自动管理生命周期。 + * + * @note 树是不可变的——节点创建后 op 和子节点关系不变，但 params 可被外部修改以支持参数优化。 + */ class SdfNode { public: - SdfOp op; - SdfParams params; - std::vector children; - std::string name; + SdfOp op; ///< 操作类型 + SdfParams params; ///< 参数值 + std::vector children; ///< 子节点 + std::string name; ///< 节点名（调试用） // ── Primitive factories ── + + /** @brief 创建球体节点 */ static SdfNodePtr sphere(double r); + + /** @brief 创建轴对齐长方体节点 */ static SdfNodePtr box(const Point3D& half_extents); + + /** @brief 创建圆角长方体节点 */ static SdfNodePtr round_box(const Point3D& half_extents, double r); + + /** @brief 创建圆柱体节点 */ static SdfNodePtr cylinder(double r, double h); + + /** @brief 创建圆环体节点 */ static SdfNodePtr torus(double major_r, double minor_r); + + /** @brief 创建胶囊体节点 */ static SdfNodePtr capsule(const Point3D& a, const Point3D& b, double r); + + /** @brief 创建圆锥台节点 */ static SdfNodePtr cone(double angle_rad, double h); + + /** @brief 创建平面节点 */ static SdfNodePtr plane(const Vector3D& normal, double offset); + + /** @brief 创建椭球体节点 */ static SdfNodePtr ellipsoid(const Point3D& radii); + + /** @brief 创建三角棱柱节点 */ static SdfNodePtr triangular_prism(double h); + + /** @brief 创建六棱柱节点 */ static SdfNodePtr hex_prism(double h); + + /** @brief 创建双环链接节点 */ static SdfNodePtr link(double r, double length, double thickness); + + /** @brief 创建楔形体节点 */ static SdfNodePtr wedge(const Point3D& extents); // ── Boolean CSG factories (two children) ── + + /** @brief 创建并集节点: A ∪ B */ static SdfNodePtr op_union(SdfNodePtr a, SdfNodePtr b); + + /** @brief 创建交集节点: A ∩ B */ static SdfNodePtr op_intersection(SdfNodePtr a, SdfNodePtr b); + + /** @brief 创建差集节点: A \ B */ static SdfNodePtr op_difference(SdfNodePtr a, SdfNodePtr b); + + /** @brief 创建平滑并集节点 */ static SdfNodePtr smooth_union(SdfNodePtr a, SdfNodePtr b, double k); + + /** @brief 创建平滑交集节点 */ static SdfNodePtr smooth_intersection(SdfNodePtr a, SdfNodePtr b, double k); + + /** @brief 创建平滑差集节点 */ static SdfNodePtr smooth_difference(SdfNodePtr a, SdfNodePtr b, double k); // ── Modifier factories (one child) ── + + /** @brief 创建圆角偏移节点 */ static SdfNodePtr round(SdfNodePtr child, double r); + + /** @brief 创建壳层修饰节点 */ static SdfNodePtr onion(SdfNodePtr child, double thickness); // ── Domain transform factories (one child) ── + + /** @brief 创建周期重复节点 */ static SdfNodePtr repeat(SdfNodePtr child, const Point3D& cell); + + /** @brief 创建 X 镜像节点 */ static SdfNodePtr mirror_x(SdfNodePtr child); + + /** @brief 创建 Y 镜像节点 */ static SdfNodePtr mirror_y(SdfNodePtr child); + + /** @brief 创建 Z 镜像节点 */ static SdfNodePtr mirror_z(SdfNodePtr child); + + /** @brief 创建平移节点 */ static SdfNodePtr translate(SdfNodePtr child, const Point3D& offset); + + /** @brief 创建旋转节点（绕 Y 轴） */ static SdfNodePtr rotate(SdfNodePtr child, double angle_rad); + + /** @brief 创建缩放节点 */ static SdfNodePtr scale(SdfNodePtr child, const Point3D& factors); + + /** @brief 创建扭曲节点 */ static SdfNodePtr twist(SdfNodePtr child, double amount); + + /** @brief 创建弯曲节点 */ static SdfNodePtr bend(SdfNodePtr child, double amount); + + /** @brief 创建拉伸节点 */ static SdfNodePtr elongate(SdfNodePtr child, const Point3D& h); + + /** @brief 创建简化弯曲节点 */ static SdfNodePtr cheap_bend(SdfNodePtr child, double amount); + + /** @brief 创建位移扰动节点 */ static SdfNodePtr displace(SdfNodePtr child, double amplitude, double frequency); - /// Visit all nodes in the tree (pre-order) + /** + * @brief 前序遍历所有节点 + * + * 对树中每个节点调用回调函数 fn，遍历顺序为前序（父节点先于子节点）。 + * + * @tparam Fn 可调用对象，签名为 void(SdfNode&) + * + * @code{.cpp} + * int leaf_count = 0; + * root->visit([&](const SdfNode& n) { + * if (n.children.empty()) ++leaf_count; + * }); + * @endcode + */ template void visit(Fn&& fn) { fn(*this); for (auto& c : children) c->visit(std::forward(fn)); } - /// Visit all nodes in the tree (const pre-order) + /** + * @brief 前序遍历（const 版本） + * + * @tparam Fn 可调用对象，签名为 void(const SdfNode&) + * @see visit(Fn&&) + */ template void visit(Fn&& fn) const { fn(*this); @@ -119,22 +282,78 @@ public: } public: + /** + * @brief 构造 SDF 节点 + * @param o 操作类型 + */ explicit SdfNode(SdfOp o) : op(o) {} }; -/// Evaluate the SDF tree at a point in world space -/// Returns signed distance: negative = inside, positive = outside +/** + * @brief 在空间点 p 处求值 SDF 树 + * + * 递归遍历表达式树，在叶节点调用对应图元的 SDF 函数， + * 在内部节点执行 CSG 操作或域变形。 + * + * **算法概要:** + * 1. 域变形节点先对 p 做变换，将变换后的点传入子节点 + * 2. CSG 布尔节点分别求值两个子节点后组合 + * 3. 图元节点直接调用内联 SDF 函数 + * + * @param root SDF 表达式树的根节点 + * @param p 世界空间中的采样点 + * @return 有符号距离（负 = 内部，零 = 表面，正 = 外部） + * + * @code{.cpp} + * auto root = SdfNode::sphere(1.0); + * double d = evaluate(root, Point3D(0.5, 0.0, 0.0)); // ≈ -0.5 + * @endcode + * + * @see estimate_bounds 估计树的包围盒 + */ [[nodiscard]] double evaluate(const SdfNodePtr& root, const Point3D& p); -/// Estimate bounding box of the SDF tree -/// Returns the half-extent from center to corner (bounding radius) +/** + * @brief 估算 SDF 树的包围半尺寸 + * + * 通过对叶节点包围盒的递归合并，估算从树根到最远角落的半径。 + * 结果是包围盒的半边长（从中心到面的距离）。 + * + * @param root SDF 树的根节点 + * @param margin 额外安全边距（默认 1.0） + * @return 包围半边长 (hx, hy, hz) + * + * @note 变形节点（旋转/弯曲/扭曲等）会使包围盒放大，以保守估计覆盖变形后区域。 + * + * @see estimate_bbox 完整 AABB 结果 + */ [[nodiscard]] Point3D estimate_bounds(const SdfNodePtr& root, double margin = 1.0); -/// Compute full AABB from estimate_bounds result +/** + * @brief SDF 树的包围盒（AABB） + * + * min/max 分别为最小和最大角点坐标，包围盒以原点对称时 min = -max。 + */ struct SdfBBox { - Point3D min; - Point3D max; + Point3D min; ///< 最小角点 + Point3D max; ///< 最大角点 }; + +/** + * @brief 估算 SDF 树的完整轴对齐包围盒 + * + * 在 estimate_bounds 基础上计算对称 AABB。 + * + * @param root SDF 树的根节点 + * @param margin 额外安全边距 + * @return AABB 包围盒 + * + * @code{.cpp} + * auto bbox = estimate_bbox(root, 0.5); + * // bbox.min = (-hx-0.5, -hy-0.5, -hz-0.5) + * // bbox.max = ( hx+0.5, hy+0.5, hz+0.5) + * @endcode + */ [[nodiscard]] SdfBBox estimate_bbox(const SdfNodePtr& root, double margin = 1.0); } // namespace vde::sdf diff --git a/include/vde/sketch/constraint_solver.h b/include/vde/sketch/constraint_solver.h index 78e0631..28935e3 100644 --- a/include/vde/sketch/constraint_solver.h +++ b/include/vde/sketch/constraint_solver.h @@ -8,72 +8,241 @@ namespace vde::sketch { using core::Point2D; using core::Vector2D; -struct SketchPoint { int id; Point2D pos; bool fixed = false; }; -struct SketchLine { int id, p0, p1; }; -struct SketchCircle { int id, center; double radius; }; +/** + * @brief 草图点元素 + * @ingroup sketch + */ +struct SketchPoint { + int id; ///< 唯一标识 + Point2D pos; ///< 2D 坐标 + bool fixed = false; ///< 是否固定（排除自由度） +}; +/** + * @brief 草图线段元素 + * @ingroup sketch + */ +struct SketchLine { + int id; ///< 唯一标识 + int p0, p1; ///< 起止点 ID +}; + +/** + * @brief 草图圆元素 + * @ingroup sketch + */ +struct SketchCircle { + int id; ///< 唯一标识 + int center; ///< 圆心点 ID + double radius; ///< 半径 +}; + +/** + * @brief 约束类型枚举 + * + * 每个约束类型对应一个残差方程，方程数见注释： + * - 距离类约束：1 个方程 + * - 位置/重合类：2 个方程（x + y） + * + * @ingroup sketch + */ enum class ConstraintType { - Horizontal, // 水平线: y1 - y0 = 0 - Vertical, // 垂直线: x1 - x0 = 0 - Parallel, // 两线平行: cross(d0, d1) = 0 - Perpendicular, // 两线垂直: dot(d0, d1) = 0 - Coincident, // 两点重合: |p1-p0| = 0 - EqualLength, // 等长度: |d0| - |d1| = 0 - Distance, // 距离约束: |p1-p0| - d = 0 - Radius, // 半径约束 - Tangent, // 相切 - Concentric, // 同心 - FixedPoint, // 固定点（硬约束，排除变量） - Fixed, // 固定点位置: x-x0=0, y-y0=0 - Angle, // 角度约束: atan2(cross(d0,d1), dot(d0,d1)) - θ = 0 + Horizontal, ///< 水平线: y1 - y0 = 0 [1 eq] + Vertical, ///< 垂直线: x1 - x0 = 0 [1 eq] + Parallel, ///< 两线平行: cross(d0, d1) = 0 [1 eq] + Perpendicular, ///< 两线垂直: dot(d0, d1) = 0 [1 eq] + Coincident, ///< 两点重合: |p1-p0| = 0 [2 eqs: Δx=0, Δy=0] + EqualLength, ///< 等长度: |d0| - |d1| = 0 [1 eq] + Distance, ///< 距离约束: |p1-p0| - d = 0 [1 eq] + Radius, ///< 半径约束 [1 eq] + Tangent, ///< 相切 [1 eq] + Concentric, ///< 同心 [2 eqs] + FixedPoint, ///< 固定点（硬约束，排除变量）[0 eq — 系统预处理] + Fixed, ///< 固定点位置: x-x0=0, y-y0=0 [2 eqs] + Angle, ///< 角度约束: atan2(cross(d0,d1), dot(d0,d1)) - θ = 0 [1 eq] }; +/** + * @brief 约束实例 + * @ingroup sketch + */ struct Constraint { - ConstraintType type; - std::vector elements; // 引用的元素索引 - double value = 0.0; // 距离/角度值 + ConstraintType type; ///< 约束类型 + std::vector elements; ///< 引用的元素索引（点/线/圆 ID） + double value = 0.0; ///< 距离/角度值 }; +/** + * @brief 求解器返回结果 + * @ingroup sketch + */ struct SolverResult { - bool converged = false; - int iterations = 0; - double residual = 0.0; - std::vector points; - std::string message; + bool converged = false; ///< 是否收敛 + int iterations = 0; ///< 实际迭代次数 + double residual = 0.0; ///< 最终残差（均方根） + std::vector points; ///< 求解后的点坐标 + std::string message; ///< 诊断信息（成功或失败原因） }; +/** + * @brief 二维几何约束求解器 + * + * 基于 Gauss-Newton 非线性最小二乘的 2D 草图约束求解器。 + * 支持点、线段、圆元素的水平/垂直/平行/垂直/相切/距离/角度等约束。 + * + * 求解流程： + * 1. 构建变量向量（非固定点的 x,y 坐标） + * 2. 对每个约束方程计算残差向量 r(x) + * 3. 数值 Jacobian J = ∂r/∂x（中心差分） + * 4. Gauss-Newton 步：Δx = -(J^T J)^{-1} J^T r + * 5. x ← x + Δx，重复直到收敛或 max_iter + * + * 自由度分析： + * DOF = 2·N_free_points - Σ(约束方程数) + * DOF > 0 → 欠约束（解不唯一） + * DOF = 0 → 恰定（应有唯一解） + * DOF < 0 → 过约束（最小二乘意义下求解） + * + * @code{.cpp} + * ConstraintSolver solver; + * int p0 = solver.add_point(0, 0); // 原点 + * int p1 = solver.add_point(100, 0); + * int p2 = solver.add_point(100, 50); + * solver.add_line(p0, p1); // 底边 + * solver.fix_point(p0); // 原点固定 + * solver.add_constraint(ConstraintType::Horizontal, {p0, p1}); + * solver.add_constraint(ConstraintType::Distance, {p0, p2}, 100.0); + * + * auto result = solver.solve(); + * if (result.converged) { + * auto p = result.points; + * } + * @endcode + * + * @ingroup sketch + */ class ConstraintSolver { public: + /** + * @brief 添加自由点 + * @param x X 坐标 + * @param y Y 坐标 + * @param fixed 是否固定（true 则排除变量） + * @return 点 ID + */ int add_point(double x, double y, bool fixed = false); + + /** + * @brief 添加线段 + * @param p0 起点 ID + * @param p1 终点 ID + * @return 线段 ID + */ int add_line(int p0, int p1); + + /** + * @brief 添加圆 + * @param center 圆心点 ID + * @param radius 半径 + * @return 圆 ID + */ int add_circle(int center, double radius); + + /** + * @brief 添加约束 + * @param type 约束类型 + * @param elements 约束涉及的元素 ID 列表 + * @param val 约束值（距离、角度等），默认 0 + * + * @note 不同类型 elements 的语义： + * Horizontal/Vertical: {p0, p1} + * Parallel/Perpendicular: {l0_p0, l0_p1, l1_p0, l1_p1} + * Distance: {p0, p1}, val = 距离 + * Angle: {l0_p0, l0_p1, l1_p0, l1_p1}, val = 角度（度） + * Radius: {circle_id}, val = 半径 + * Tangent: {line_p0, line_p1, circle_id} + * Coincident: {p0, p1} + */ void add_constraint(ConstraintType type, const std::vector& elements, double val = 0.0); + + /** + * @brief 标记点为固定（硬约束） + * @param pid 点 ID + * @note 固定点不参与变量优化，从自由度中排除 + */ void fix_point(int pid) { if (pid >= 0 && pid < static_cast(points_.size())) points_[pid].fixed = true; } + /** + * @brief 执行约束求解 + * @param max_iter 最大 Gauss-Newton 迭代次数，默认 100 + * @param tol 收敛容差（梯度范数 < tol 时停止），默认 1e-8 + * @return SolverResult 求解结果 + * + * @note 收敛判断：||J^T r|| < tol + * @note 欠约束时（DOF > 0）解不唯一；过约束时（DOF < 0）为最小二乘解 + */ [[nodiscard]] SolverResult solve(int max_iter = 100, double tol = 1e-8); + + /** + * @brief 计算自由度 + * @return DOF = 2·N_free - Σ(约束方程数) + * + * @note DOF < 0 时系统过约束，求解器仍然工作但解非精确满足所有约束 + */ [[nodiscard]] int degrees_of_freedom() const; + + /** + * @brief 当前点集访问 + * @return SketchPoint 数组只读引用 + */ [[nodiscard]] const std::vector& points() const { return points_; } + + /** + * @brief 查询点坐标 + * @param id 点 ID + * @return Point2D 当前坐标 + */ [[nodiscard]] Point2D get_point(int id) const; private: - std::vector points_; - std::vector lines_; - std::vector circles_; - std::vector constraints_; + std::vector points_; ///< 点集 + std::vector lines_; ///< 线段集 + std::vector circles_; ///< 圆集 + std::vector constraints_; ///< 约束集 - // 辅助：同步 vars → points_（用于 Jacobian 数值计算） + /** + * @brief 将求解器变量向量同步回 points_（用于 Jacobian 数值计算） + * @param vars 变量数组（长度为 2·N_free） + * @param nv 变量数 + */ void sync_from_vars(const double* vars, int nv); - // 约束残差评估 + /** + * @brief 评估单个约束的残差值 + * @param c 约束 + * @return 残差（方程 → 0 表示约束满足） + */ double eval_constraint(const Constraint& c) const; - // 约束方程数（Coincident/Fixed → 2，其他 → 1） + /** + * @brief 计算单个约束的方程数 + * @param c 约束 + * @return Coincident/Fixed → 2，其余 → 1 + */ int count_equations(const Constraint& c) const; - // Jacobian 填充：将约束 ci 的 Jacobian 行写入矩阵 J + /** + * @brief 填充 Jacobian 矩阵中对应约束 ci 的行 + * @param c 约束 + * @param ci 约束索引 + * @param row_map 方程行映射 + * @param J_data Jacobian 矩阵数据（行优先） + * @param nv 变量数 + * @param stride 行跨度 + */ void fill_jacobian_rows(const Constraint& c, int ci, const std::vector& row_map, double* J_data, int nv, int stride) const;