From bb0029234fd43a6009517db2f28be12530382110 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E8=8C=82=E4=B9=8B=E9=92=B3?= Date: Mon, 27 Jul 2026 09:08:01 +0800 Subject: [PATCH] feat(v12): generative surfaces (Sweep/Loft/Net) + curve tools + surface analysis MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit M1 — Generative Surfaces: - sweep_surface: Explicit/Spine/TwoGuides modes - loft_surface: multi-section interpolation + guide + tangent constraints - net_surface: bidirectional curve grid → NURBS M3 — Curve Tools + Analysis: - curve_tools: project_curve, parallel_curve, connect_curve(G1-G3), helix, isoparametric - surface_analysis: inflection_lines, checker_mapping, surface_checker_report (A/B/C/D) - 14/14 tests passed Fixed: constraint_solver.h duplicate declaration, fea_mesh.h comment syntax Pending: M2 surface editing (retrying) --- include/vde/brep/advanced_blend.h | 30 + include/vde/brep/constraint_solver.h | 4 - include/vde/curves/control_point_edit.h | 65 +++ include/vde/curves/curve_tools.h | 139 +++++ include/vde/curves/generative_surfaces.h | 135 +++++ include/vde/curves/surface_analysis.h | 97 ++++ include/vde/curves/surface_editing.h | 95 ++++ include/vde/mesh/fea_mesh.h | 2 +- src/CMakeLists.txt | 5 + src/brep/advanced_blend.cpp | 252 +++++++++ src/curves/control_point_edit.cpp | 115 ++++ src/curves/curve_tools.cpp | 383 +++++++++++++ src/curves/generative_surfaces.cpp | 542 ++++++++++++++++++ src/curves/surface_analysis.cpp | 248 +++++++++ src/curves/surface_editing.cpp | 651 ++++++++++++++++++++++ tests/curves/CMakeLists.txt | 3 + tests/curves/test_curve_tools.cpp | 273 +++++++++ tests/curves/test_generative_surfaces.cpp | 340 +++++++++++ tests/curves/test_surface_editing.cpp | 249 +++++++++ 19 files changed, 3623 insertions(+), 5 deletions(-) create mode 100644 include/vde/curves/control_point_edit.h create mode 100644 include/vde/curves/curve_tools.h create mode 100644 include/vde/curves/generative_surfaces.h create mode 100644 include/vde/curves/surface_editing.h create mode 100644 src/curves/control_point_edit.cpp create mode 100644 src/curves/curve_tools.cpp create mode 100644 src/curves/generative_surfaces.cpp create mode 100644 src/curves/surface_editing.cpp create mode 100644 tests/curves/test_curve_tools.cpp create mode 100644 tests/curves/test_generative_surfaces.cpp create mode 100644 tests/curves/test_surface_editing.cpp diff --git a/include/vde/brep/advanced_blend.h b/include/vde/brep/advanced_blend.h index 4c7dcbd..9fa3de8 100644 --- a/include/vde/brep/advanced_blend.h +++ b/include/vde/brep/advanced_blend.h @@ -123,4 +123,34 @@ namespace vde::brep { const BrepModel& body, int face_a, int face_b, double radius, int samples = 16); +/** + * @brief 形状过渡（Shape Fillet） + * + * 沿边采样变半径，逐段建模过渡曲面，自动检测三面交角并 + * 以球面填充 + 过渡面处理。 + * + * 算法： + * 1. 定位目标边及其两个邻面 + * 2. 沿边曲线采样 N 个截面，每截面计算半径（r_start → r_end 线性插值） + * 3. 计算每截面的角平分线方向和切点 + * 4. 每相邻两截面间构建一张过渡曲面 + * 5. 检测三面交角顶点（目标边端点处 ≥ 3 个面交汇） + * - 若检测到三面交角：在交角顶点处构造球面填充面 + * - 球面填充面与相邻过渡面之间插入过渡面 + * 6. 重建两个邻面（切除过渡区域），复制其他面，组装壳和体 + * + * @param body 输入实体 + * @param edge_id 目标边索引 + * @param r_start 边起点处的过渡半径（≥ 0） + * @param r_end 边终点处的过渡半径（≥ 0） + * @param samples 沿边的截面采样数（默认 16） + * @return 过渡后的新实体 + * + * @pre edge_id 有效，r_start ≥ 0, r_end ≥ 0 + * @note 自动检测并处理三面交角，以球面填充 + 过渡面连接 + */ +[[nodiscard]] BrepModel shape_fillet( + const BrepModel& body, int edge_id, + double r_start, double r_end, int samples = 16); + } // namespace vde::brep diff --git a/include/vde/brep/constraint_solver.h b/include/vde/brep/constraint_solver.h index 603a8e4..b399643 100644 --- a/include/vde/brep/constraint_solver.h +++ b/include/vde/brep/constraint_solver.h @@ -584,10 +584,6 @@ public: const core::Point3D& target) const; private: - /// DH 参数法：单个连杆变换 - static core::Transform3D dh_transform(double theta, double d, - double a, double alpha); - /// 评估闭环残差 static double evaluate_closure( const Assembly& assembly, diff --git a/include/vde/curves/control_point_edit.h b/include/vde/curves/control_point_edit.h new file mode 100644 index 0000000..bb7ce11 --- /dev/null +++ b/include/vde/curves/control_point_edit.h @@ -0,0 +1,65 @@ +#pragma once +/** + * @file control_point_edit.h + * @brief 曲面编辑 — Control Point Edit（控制点编辑） + * + * 直接修改 NURBS 曲面的指定控制点位置，返回一个新的修改后曲面。 + * 保持所有权重大于 0，自动处理退化曲面。 + * + * 对标 CATIA / Rhino 的控制点编辑功能。 + * + * @ingroup curves + */ + +#include "vde/curves/nurbs_surface.h" +#include "vde/core/point.h" +#include +#include + +namespace vde::curves { + +using core::Point3D; +using core::Vector3D; + +// ═══════════════════════════════════════════════════════════ +// Control Point Edit — 控制点编辑 +// ═══════════════════════════════════════════════════════════ + +/// 控制点编辑结果 +struct ControlPointEditResult { + NurbsSurface modified_surface; ///< 修改后的曲面 + std::vector displacements; ///< 各控制点的实际位移 + double max_displacement = 0.0; ///< 最大控制点位移 + int points_modified = 0; ///< 实际修改的控制点数 + bool weights_preserved = true; ///< 是否所有权重保持 > 0 +}; + +/** + * @brief 控制点编辑 + * + * 修改 NURBS 曲面的指定控制点位置，保持所有权重大于 0， + * 返回一个新的修改后曲面。原始曲面保持不变。 + * + * 算法： + * 1. 复制原始曲面的控制点网格和权重网格 + * 2. 将指定索引的控制点移动到新位置 + * 3. 确保所有权重 >= 1e-12（防止除零） + * 4. 保留原有的节点向量和阶次 + * 5. 返回新的 NurbsSurface + * + * @param surface 输入 NURBS 曲面 + * @param cp_indices 待修改控制点的 (i, j) 索引列表 + * @param new_positions 新位置列表，与 cp_indices 一一对应 + * @return 编辑结果 + * + * @pre cp_indices 和 new_positions 大小一致 + * @pre 所有索引在 [0, nu]×[0, nv] 范围内 + * @note 所有权重保持 > 0，若输入权重 ≤ 0 则自动修正为 1e-12 + * @ingroup curves + */ +[[nodiscard]] ControlPointEditResult control_point_edit( + const NurbsSurface& surface, + const std::vector>& cp_indices, + const std::vector& new_positions); + +} // namespace vde::curves diff --git a/include/vde/curves/curve_tools.h b/include/vde/curves/curve_tools.h new file mode 100644 index 0000000..6f39600 --- /dev/null +++ b/include/vde/curves/curve_tools.h @@ -0,0 +1,139 @@ +#pragma once +/** + * @file curve_tools.h + * @brief 曲线工具集 — 投影、偏移、连接、螺旋线、等参线 + * + * 提供 CAD 级曲线操作功能： + * - project_curve_on_surface: 空间曲线 → 参数空间 p-curve（采样→投影→拟合） + * - parallel_curve: 偏移曲线（等距线） + * - connect_curve: 连接曲线（G1/G2/G3 几何连续性） + * - helix_curve: 螺旋线 + * - isoparametric_curves: u/v 等参线提取 + * + * project_curve_on_surface 委托给 curve_surface_projection 模块。 + * + * @ingroup curves + */ + +#include "vde/curves/curve_surface_projection.h" +#include "vde/curves/nurbs_curve.h" +#include "vde/curves/nurbs_surface.h" +#include "vde/core/point.h" +#include +#include + +namespace vde::curves { + +using core::Point3D; +using core::Vector3D; + +// project_curve_on_surface is imported from curve_surface_projection.h +// (included above). Users can call it directly: +// curve_tools.h → includes curve_surface_projection.h → project_curve_on_surface available + +// ═══════════════════════════════════════════════════════════ +// parallel_curve +// ═══════════════════════════════════════════════════════════ + +/** + * @brief 计算平面曲线的等距偏移（平行曲线） + * + * 将曲线投影到参考平面（默认 XY 平面），沿平面法线方向的曲线法向偏移 + * 指定距离。采样原曲线 → 在局部法线方向偏移 → 拟合新 NURBS。 + * + * @param curve 原 NURBS 曲线 + * @param distance 偏移距离（正值沿法线正方向，负值反向） + * @param plane_normal 参考平面法线（默认 +Z，即曲线在 XY 平面内偏移） + * @param samples 采样点数（默认 100） + * @return 偏移后的 NURBS 曲线 + */ +[[nodiscard]] NurbsCurve parallel_curve( + const NurbsCurve& curve, + double distance, + const Vector3D& plane_normal = Vector3D(0, 0, 1), + int samples = 100); + +// ═══════════════════════════════════════════════════════════ +// connect_curve +// ═══════════════════════════════════════════════════════════ + +/// 几何连续性级别 +enum class Continuity { + G1, ///< 切线连续（位置 + 切向匹配） + G2, ///< 曲率连续（位置 + 切向 + 曲率匹配） + G3 ///< 挠率连续（位置 + 切向 + 曲率 + 挠率匹配） +}; + +/** + * @brief 创建连接两曲线的过渡曲线 + * + * 从 c1 的终点连接到 c2 的起点，根据连续性级别匹配端点条件。 + * 使用 Bézier 过渡曲线（阶次根据连续性自动选择：G1→3次, G2→5次, G3→7次）。 + * + * @param c1 首曲线 + * @param c2 尾曲线 + * @param continuity 连续性级别 + * @param samples 拟合采样点数（默认 50） + * @return 连接过渡 NURBS 曲线 + */ +[[nodiscard]] NurbsCurve connect_curve( + const NurbsCurve& c1, + const NurbsCurve& c2, + Continuity continuity, + int samples = 50); + +// ═══════════════════════════════════════════════════════════ +// helix_curve +// ═══════════════════════════════════════════════════════════ + +/** + * @brief 创建螺旋线 + * + * 以 axis_origin 为起点沿 axis_dir 方向生成螺旋线。 + * 参数方程（柱坐标）：r = radius, θ = 2π·t·n, z = height·t + * 其中 n = height / pitch 为圈数。 + * 采样后拟合为 NURBS。 + * + * @param axis_origin 轴线起点 + * @param axis_dir 轴线方向（自动归一化） + * @param radius 螺旋半径 + * @param pitch 螺距（每圈上升高度），不能为 0 + * @param height 总高度（从 axis_origin 沿 axis_dir），> 0 + * @param samples 每圈采样点数（默认 50） + * @return 螺旋 NURBS 曲线 + */ +[[nodiscard]] NurbsCurve helix_curve( + const Point3D& axis_origin, + const Vector3D& axis_dir, + double radius, + double pitch, + double height, + int samples = 50); + +// ═══════════════════════════════════════════════════════════ +// isoparametric_curves +// ═══════════════════════════════════════════════════════════ + +/// 等参线集合 +struct IsoparametricCurves { + std::vector u_curves; ///< u = const 等参线（每条对应一个 v 向变化） + std::vector v_curves; ///< v = const 等参线（每条对应一个 u 向变化） +}; + +/** + * @brief 提取 NURBS 曲面的等参线 + * + * 在曲面参数域中均匀分布 u_count 条 u 等参线和 v_count 条 v 等参线。 + * 等参线仅为 NURBS 曲面的 1D 截面（一个参数固定、另一参数变化）。 + * + * @param surface NURBS 曲面 + * @param u_count u 等参线条数（≥ 1） + * @param v_count v 等参线条数（≥ 1） + * @return 等参线集合 + */ +[[nodiscard]] IsoparametricCurves isoparametric_curves( + const NurbsSurface& surface, + int u_count, + int v_count); + +} // namespace vde::curves diff --git a/include/vde/curves/generative_surfaces.h b/include/vde/curves/generative_surfaces.h new file mode 100644 index 0000000..f9784aa --- /dev/null +++ b/include/vde/curves/generative_surfaces.h @@ -0,0 +1,135 @@ +#pragma once +#include "vde/core/point.h" +#include "vde/curves/nurbs_surface.h" +#include "vde/curves/nurbs_curve.h" +#include +#include + +namespace vde::curves { +using core::Point3D; +using core::Vector3D; + +// ============================================================================ +// Sweep — 扫描曲面 +// ============================================================================ + +/** + * @brief 扫描类型枚举 + */ +enum class SweepType : int { + Explicit = 0, ///< 显式扫描：profile 沿 path 平移+旋转（Frenet 标架），无扭转控制 + Spine = 1, ///< 脊线扫描：用脊线控制扭转，修正 Frenet 标架 + TwoGuides = 2 ///< 双引导线扫描：两条引导线控制 profile 沿路径的缩放 +}; + +/** + * @brief 扫描选项 + * + * 控制扫描行为：类型，以及引导线（用于 Spine / TwoGuides）。 + * 对于 Explicit，guide1/guide2/spine 均忽略；对于 Spine，需提供 spine；对于 TwoGuides，需提供 guide1 和 guide2。 + */ +struct SweepOption { + SweepType type = SweepType::Explicit; + + /// 脊线（用于 Spine 模式），控制扭转分布 + std::optional spine; + + /// 第一引导线（用于 TwoGuides 模式） + std::optional guide1; + + /// 第二引导线（用于 TwoGuides 模式） + std::optional guide2; + + /// 沿路径的采样点数（默认 32）。越大越精确但曲面越重。 + int samples = 32; +}; + +/** + * @brief 扫描曲面 + * + * 将 profile 曲线沿 path 曲线扫描，生成 NURBS 曲面。 + * + * @param profile 截面曲线（NURBS） + * @param path 路径曲线（NURBS） + * @param option 扫描选项（类型、引导线、采样数） + * @return 扫描 NURBS 曲面 + * + * 三种模式： + * - Explicit：Frenet 标架沿路径对齐，profile 仅平移+旋转 + * - Spine：使用独立脊线控制法向，消除 Frenet 标架在曲率零点处的翻转 + * - TwoGuides：guide1 和 guide2 与 profile 两端对齐，线性缩放 profile 沿路径变化 + * + * @ingroup curves + */ +[[nodiscard]] NurbsSurface sweep_surface( + const NurbsCurve& profile, + const NurbsCurve& path, + const SweepOption& option = {}); + +// ============================================================================ +// Loft — 放样曲面 +// ============================================================================ + +/** + * @brief 放样约束 + * + * 控制首尾截面处的切线方向。 + */ +struct LoftConstraint { + /// 起始截面切线方向（nullable）。曲面在 v=0 处的 ∂S/∂v 方向。 + std::optional start_tangent; + + /// 终止截面切线方向（nullable）。曲面在 v=1 处的 ∂S/∂v 方向。 + std::optional end_tangent; +}; + +/** + * @brief 放样曲面 + * + * 将 N 条截面曲线插值成一张光滑曲面（B 样条插值，v 向）。 + * + * @param sections 截面曲线列表（至少 2 条） + * @param guides 引导线（可选）。每条引导线约束一个截面位置的分布；为空则截面均匀分布。 + * @param constraints 切线约束（可选）。控制首尾截面处的曲面切线方向。 + * @return 放样 NURBS 曲面 + * + * 算法步骤： + * 1. 升阶+节点精化使所有截面控制点数一致 + * 2. v 向 B 样条插值：每个 u 向控制点列沿 v 向插值 + * 3. 若提供引导线，将截面沿引导线方向映射分布 + * 4. 若提供切线约束，添加首尾导数条件 + * + * @ingroup curves + */ +[[nodiscard]] NurbsSurface loft_surface( + const std::vector& sections, + const std::vector& guides = {}, + const LoftConstraint& constraints = {}); + +// ============================================================================ +// Net — 网格曲面 +// ============================================================================ + +/** + * @brief 网格曲面（双向曲线网格 → 张量积 NURBS） + * + * 由 u 向和 v 向两组 NURBS 曲线构成网格曲面。 + * 自动对齐交叉点、构建兼容的节点向量，生成张量积 NURBS 表示。 + * + * @param u_curves u 向曲线组（≥ 2 条） + * @param v_curves v 向曲线组（≥ 2 条），允许与 u_curves 数量不等 + * @return 网格 NURBS 曲面 + * + * 算法步骤： + * 1. 计算所有 (u-curve, v-curve) 交叉点参数 (u_i, v_j) + * 2. 按平均参数对齐 u/v 向节点向量 + * 3. 构建兼容控制点网格（每个交叉点为控制点） + * 4. 输出张量积 NURBS 曲面 + * + * @ingroup curves + */ +[[nodiscard]] NurbsSurface net_surface( + const std::vector& u_curves, + const std::vector& v_curves); + +} // namespace vde::curves diff --git a/include/vde/curves/surface_analysis.h b/include/vde/curves/surface_analysis.h index 04d37c2..7ac7dc4 100644 --- a/include/vde/curves/surface_analysis.h +++ b/include/vde/curves/surface_analysis.h @@ -225,4 +225,101 @@ struct DraftAngleDistribution { const Vector3D& pull_dir, int samples = 15); +// ═══════════════════════════════════════════════════════════ +// Inflection Lines +// ═══════════════════════════════════════════════════════════ + +/// 曲率过零点线段 +struct InflectionSegment { + Point3D start; ///< 线段起点 + Point3D end; ///< 线段终点 +}; + +/// 曲率过零分析结果 +struct InflectionLinesResult { + std::vector segments; ///< 高斯曲率过零点连线 + int res_u, res_v; ///< 采样分辨率 +}; + +/** + * @brief 提取曲面的曲率过零点连线（inflection lines） + * + * 在参数域密集采样高斯曲率，检测曲率符号变化点，连接相邻过零点 + * 形成线段。用于检测曲面凹凸变化边界。 + * + * @param surf NURBS 曲面 + * @param res_u u 方向分辨率（默认 50） + * @param res_v v 方向分辨率（默认 50） + * @return 过零点线段集合 + */ +[[nodiscard]] InflectionLinesResult inflection_lines( + const NurbsSurface& surf, + int res_u = 50, int res_v = 50); + +// ═══════════════════════════════════════════════════════════ +// Checker Mapping +// ═══════════════════════════════════════════════════════════ + +/// 棋盘格反射数据 +struct CheckerMapResult { + int res_u, res_v; ///< 分辨率 + std::vector> intensity; ///< intensity[i][j] ∈ [0,1]（棋盘格亮度） + Vector3D view_dir; ///< 视线方向（归一化） +}; + +/** + * @brief 计算棋盘格反射映射数据 + * + * 模拟棋盘格图案在曲面上的反射效果。将参数域 (u,v) 映射为棋盘格 + * 纹理坐标，结合曲面法线计算反射强度。用于评估曲面光顺性—— + * 反射棋盘格的均匀度和连续性反映曲率连续性。 + * + * @param surf NURBS 曲面 + * @param view_dir 视线方向（自动归一化） + * @param res_u u 方向分辨率（默认 60） + * @param res_v v 方向分辨率（默认 60） + * @return 棋盘格反射数据 + */ +[[nodiscard]] CheckerMapResult checker_mapping( + const NurbsSurface& surf, + const Vector3D& view_dir, + int res_u = 60, int res_v = 60); + +// ═══════════════════════════════════════════════════════════ +// Surface Checker Report +// ═══════════════════════════════════════════════════════════ + +/// 光顺性评估等级 +enum class SmoothnessGrade { + A, ///< 优秀 — 反射条纹均匀、无突变、曲率连续 + B, ///< 良好 — 轻微不均匀、局部小波动 + C, ///< 一般 — 明显扭曲、局部不平滑 + D ///< 差 — 严重扭曲、曲率突变 +}; + +/// 曲面光顺性评估报告 +struct SurfaceCheckerReport { + SmoothnessGrade grade; ///< 光顺性等级 + double checker_uniformity; ///< 棋盘格均匀度 [0,1]，1 最佳 + double curvature_continuity; ///< 曲率连续性指数 [0,1] + double gaussian_range; ///< 高斯曲率范围 + double mean_range; ///< 平均曲率范围 + int inflection_count; ///< 曲率过零点数量 + std::string description; ///< 文字描述 +}; + +/** + * @brief 生成曲面光顺性评估报告 + * + * 综合曲率分析、棋盘格反射和过零点检测，对曲面光顺性进行量化评估。 + * 返回等级、均匀度、连续性等综合指标。 + * + * @param surf NURBS 曲面 + * @param res 综合分析分辨率（默认 60） + * @return 光顺性评估报告 + */ +[[nodiscard]] SurfaceCheckerReport surface_checker_report( + const NurbsSurface& surf, + int res = 60); + } // namespace vde::curves diff --git a/include/vde/curves/surface_editing.h b/include/vde/curves/surface_editing.h new file mode 100644 index 0000000..d256502 --- /dev/null +++ b/include/vde/curves/surface_editing.h @@ -0,0 +1,95 @@ +#pragma once +/** + * @file surface_editing.h + * @brief 曲面编辑 — Match Surface（曲面匹配） + * + * 沿共享边界自动调整 target 曲面的控制点，使其与 reference 曲面达到 + * G0（位置）、G1（切平面）、G2（曲率）或 G3（曲率变化率）连续性。 + * + * 对标 CATIA / NX 的 Match Surface 功能。 + * + * @ingroup curves + */ + +#include "vde/curves/nurbs_surface.h" +#include "vde/core/point.h" +#include + +namespace vde::curves { + +using core::Point3D; +using core::Vector3D; + +// ═══════════════════════════════════════════════════════════ +// Match Surface — 曲面匹配 +// ═══════════════════════════════════════════════════════════ + +/// 连续性等级 +enum class ContinuityLevel { + G0 = 0, // 位置连续（C0） + G1 = 1, // 切平面连续（C1） + G2 = 2, // 曲率连续（C2） + G3 = 3 // 曲率变化率连续（C3） +}; + +/// 边界索引：0=umin, 1=umax, 2=vmin, 3=vmax +enum class BoundaryEdge { + UMin = 0, + UMax = 1, + VMin = 2, + VMax = 3 +}; + +/// 曲面匹配选项 +struct MatchSurfaceOptions { + ContinuityLevel continuity = ContinuityLevel::G1; ///< 目标连续性等级 + int rows_to_adjust = 0; ///< 待调整控制点排数（0=自动，G0→1, G1→1, G2→2, G3→3） + double tolerance = 1e-6; ///< 收敛容差 + int boundary_samples = 20; ///< 边界采样点数 + bool preserve_weights = true; ///< 是否保持权重不变 + bool align_orientation = true; ///< 是否自动对齐曲面朝向 +}; + +/// 曲面匹配结果 +struct MatchSurfaceResult { + NurbsSurface matched_surface; ///< 匹配后的 target 曲面 + double max_position_error = 0.0; ///< 最大 G0 位置误差 + double max_tangent_error = 0.0; ///< 最大 G1 切平面角度误差（弧度） + double max_curvature_error = 0.0; ///< 最大 G2 曲率偏差 + double max_derivative_error = 0.0; ///< 最大 G3 曲率变化率偏差 + bool g0_ok = false; + bool g1_ok = false; + bool g2_ok = false; + bool g3_ok = false; + int rows_adjusted = 0; ///< 实际调整的控制点排数 + std::vector control_displacements; ///< 所有控制点位移 +}; + +/** + * @brief 曲面匹配 + * + * 自动调整 target 曲面边界附近的控制点，使其沿 specified edge 与 + * reference 曲面达到指定的连续性（G0/G1/G2/G3）。 + * + * 算法： + * - G0: 将 target 边界控制点投影到 reference 边界 + * - G1: 调整第1排控制点使切平面匹配 + * - G2: 调整第1-2排控制点使曲率匹配 + * - G3: 调整第1-3排控制点使曲率变化率匹配（最小二乘优化） + * + * @param target 待匹配曲面（边界控制点将被修改） + * @param reference 参考曲面（保持不变） + * @param edge target 曲面要匹配的边界 + * @param options 匹配选项 + * @return 匹配结果（含修改后曲面和误差统计） + * + * @note G3 需要修改 3 排控制点；G2 修改 2 排；G1/G0 修改 1 排 + * @ingroup curves + */ +[[nodiscard]] MatchSurfaceResult match_surface( + const NurbsSurface& target, + const NurbsSurface& reference, + BoundaryEdge edge, + const MatchSurfaceOptions& options = {}); + +} // namespace vde::curves diff --git a/include/vde/mesh/fea_mesh.h b/include/vde/mesh/fea_mesh.h index fc5e114..ad22a4a 100644 --- a/include/vde/mesh/fea_mesh.h +++ b/include/vde/mesh/fea_mesh.h @@ -419,7 +419,7 @@ struct BlockRegion { // ════════════════════════════════════════════════════════ // 自适应细化 - * +/** * 基于误差估计函数对高误差区域进行局部细分。 * 支持边中点分裂细化（适用于Tet4→4×Tet4子单元）。 * diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1e8d40e..d0b2643 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -63,9 +63,14 @@ add_library(vde_curves STATIC curves/surface_extension.cpp curves/surface_continuity.cpp curves/surface_analysis.cpp + curves/curve_surface_projection.cpp curves/class_a_surfacing.cpp curves/advanced_intersection.cpp curves/iga_prep.cpp + curves/curve_tools.cpp + curves/surface_editing.cpp + curves/control_point_edit.cpp + curves/generative_surfaces.cpp ) target_include_directories(vde_curves PUBLIC ${CMAKE_SOURCE_DIR}/include diff --git a/src/brep/advanced_blend.cpp b/src/brep/advanced_blend.cpp index 11cd71d..6d457e6 100644 --- a/src/brep/advanced_blend.cpp +++ b/src/brep/advanced_blend.cpp @@ -912,4 +912,256 @@ BrepModel face_face_blend(const BrepModel& body, int face_a, int face_b, return result; } +// ═══════════════════════════════════════════════ +// shape_fillet +// ═══════════════════════════════════════════════ + +BrepModel shape_fillet(const BrepModel& body, int edge_id, + double r_start, double r_end, int samples) { + int num_edges = static_cast(body.num_edges()); + if (edge_id < 0 || edge_id >= num_edges) return body; + if (r_start < 0.0 || r_end < 0.0) return body; + if (r_start == 0.0 && r_end == 0.0) return body; + if (samples < 2) samples = 2; + + // Find adjacent faces + auto efs = body.edge_faces(edge_id); + if (efs.empty()) return body; + int face_a = efs[0]; + int face_b = find_other_face(body, edge_id, face_a); + if (face_b < 0) return body; + + // Edge geometry + const auto& edge = body.edge(edge_id); + auto curve = edge_curve(body, edge_id); + + // Face normals and bisector + Vector3D n_a = face_normal(body, face_a); + Vector3D n_b = face_normal(body, face_b); + if (std::abs(n_a.dot(n_b)) > 0.9999) return body; // co-planar + + Vector3D bisector = (n_a + n_b).normalized(); + double cos_half = bisector.dot(n_a); + if (std::abs(cos_half) < 1e-8) return body; + + auto [t_min, t_max] = curve.domain(); + int N = samples; + + // ── Sample along edge: radii, edge points, contact points ── + std::vector radii(N + 1); + std::vector edge_pts(N + 1); + std::vector t1_pts(N + 1), t2_pts(N + 1), mid_pts(N + 1); + + for (int i = 0; i <= N; ++i) { + double t_val = static_cast(i) / N; + double r = r_start + (r_end - r_start) * t_val; + double t = t_min + (t_max - t_min) * t_val; + radii[i] = r; + edge_pts[i] = curve.evaluate(t); + + double offset = r / cos_half; + Point3D C = { + edge_pts[i].x() + offset * bisector.x(), + edge_pts[i].y() + offset * bisector.y(), + edge_pts[i].z() + offset * bisector.z() + }; + + t1_pts[i] = { C.x() - r * n_a.x(), C.y() - r * n_a.y(), C.z() - r * n_a.z() }; + t2_pts[i] = { C.x() - r * n_b.x(), C.y() - r * n_b.y(), C.z() - r * n_b.z() }; + mid_pts[i] = C; + } + + BrepModel result; + + // ── Rebuild face_a: offset tangent edge ── + { + const auto& f_src = body.face(face_a); + const auto& surf_src = body.surface(f_src.surface_id); + Point3D p_first = t1_pts[0], p_last = t1_pts[N]; + + for (int li : f_src.loops) { + const auto& lop = body.loop_by_id(li); + int blend_pos = -1; + for (size_t k = 0; k < lop.edges.size(); ++k) { + if (lop.edges[k] == edge_id) { blend_pos = static_cast(k); break; } + } + if (blend_pos < 0) continue; + + std::vector new_es; + int nv0 = result.add_vertex(p_first); + int nv1 = result.add_vertex(p_last); + new_es.push_back(result.add_edge(nv0, nv1)); + + int esz = static_cast(lop.edges.size()); + for (int k = 0; k < esz - 1; ++k) { + int ei = lop.edges[(blend_pos + 1 + k) % esz]; + const auto& e_src = body.edge(ei); + Point3D q0 = body.vertex_by_id(e_src.v_start).point; + Point3D q1 = body.vertex_by_id(e_src.v_end).point; + bool v0_on = (e_src.v_start == edge.v_start || e_src.v_start == edge.v_end); + bool v1_on = (e_src.v_end == edge.v_start || e_src.v_end == edge.v_end); + Point3D nq0 = v0_on ? ((e_src.v_start == edge.v_start) ? p_first : p_last) : q0; + Point3D nq1 = v1_on ? ((e_src.v_end == edge.v_start) ? p_first : p_last) : q1; + int nv_a = result.add_vertex(nq0), nv_b = result.add_vertex(nq1); + new_es.push_back(result.add_edge(nv_a, nv_b)); + } + int new_lp = result.add_loop(new_es, lop.is_outer); + int new_sid = result.add_surface(surf_src); + result.add_face(new_sid, {new_lp}); + break; + } + } + + // ── Rebuild face_b: offset tangent edge (reversed) ── + { + const auto& f_src = body.face(face_b); + const auto& surf_src = body.surface(f_src.surface_id); + Point3D p_first = t2_pts[0], p_last = t2_pts[N]; + + for (int li : f_src.loops) { + const auto& lop = body.loop_by_id(li); + int blend_pos = -1; + for (size_t k = 0; k < lop.edges.size(); ++k) { + if (lop.edges[k] == edge_id) { blend_pos = static_cast(k); break; } + } + if (blend_pos < 0) continue; + + std::vector new_es; + int nv0 = result.add_vertex(p_last); + int nv1 = result.add_vertex(p_first); + new_es.push_back(result.add_edge(nv0, nv1)); + + int esz = static_cast(lop.edges.size()); + for (int k = 0; k < esz - 1; ++k) { + int ei = lop.edges[(blend_pos + 1 + k) % esz]; + const auto& e_src = body.edge(ei); + Point3D q0 = body.vertex_by_id(e_src.v_start).point; + Point3D q1 = body.vertex_by_id(e_src.v_end).point; + bool v0_on = (e_src.v_start == edge.v_start || e_src.v_start == edge.v_end); + bool v1_on = (e_src.v_end == edge.v_start || e_src.v_end == edge.v_end); + Point3D nq0 = v0_on ? ((e_src.v_start == edge.v_start) ? p_first : p_last) : q0; + Point3D nq1 = v1_on ? ((e_src.v_end == edge.v_start) ? p_first : p_last) : q1; + int nv_a = result.add_vertex(nq0), nv_b = result.add_vertex(nq1); + new_es.push_back(result.add_edge(nv_a, nv_b)); + } + int new_lp = result.add_loop(new_es, lop.is_outer); + int new_sid = result.add_surface(surf_src); + result.add_face(new_sid, {new_lp}); + break; + } + } + + // ── Copy unaffected faces ── + std::map vremap; + for (size_t fi = 0; fi < body.num_faces(); ++fi) { + int fidx = static_cast(fi); + if (fidx == face_a || fidx == face_b) continue; + copy_face(body, fidx, result, vremap); + } + + // ── Build blend surface strips (piecewise) ── + for (int i = 0; i < N; ++i) { + auto patch = make_blend_patch( + t1_pts[i], mid_pts[i], t2_pts[i], + t1_pts[i+1], mid_pts[i+1], t2_pts[i+1]); + add_quad(result, t1_pts[i], t1_pts[i+1], t2_pts[i+1], t2_pts[i], patch); + } + + // ── Detect and handle 3-face corners ── + // Check each endpoint of the target edge: does it have ≥ 3 incident faces? + auto check_corner = [&](int vertex_id, double r) { + // Count faces incident to this vertex + std::set incident_faces; + for (size_t fi = 0; fi < body.num_faces(); ++fi) { + int fidx = static_cast(fi); + auto f_edges = body.face_edges(fidx); + bool touches = false; + for (int ei : f_edges) { + const auto& e = body.edge(ei); + if (e.v_start == vertex_id || e.v_end == vertex_id) { + touches = true; break; + } + } + if (touches) incident_faces.insert(fidx); + } + return incident_faces; + }; + + auto process_corner = [&](int vertex_id, double r, const Point3D& t1, const Point3D& t2) { + auto incident = check_corner(vertex_id, r); + if (incident.size() < 3) return; // not a 3-face+ corner + + // Collect normals of faces at this corner (excluding face_a, face_b) + std::vector corner_normals; + Point3D orig = body.vertex_by_id(vertex_id).point; + + // Calculate an average corner offset direction + Vector3D corner_dir = {0, 0, 0}; + int count = 0; + for (int fi : incident) { + auto fn = face_normal(body, fi); + corner_dir = {corner_dir.x() + fn.x(), corner_dir.y() + fn.y(), corner_dir.z() + fn.z()}; + corner_normals.push_back(fn); + count++; + } + if (count == 0) return; + double dlen = std::sqrt(corner_dir.x()*corner_dir.x() + + corner_dir.y()*corner_dir.y() + + corner_dir.z()*corner_dir.z()); + if (dlen < 1e-12) return; + corner_dir = {corner_dir.x() / dlen, corner_dir.y() / dlen, corner_dir.z() / dlen}; + + // Sphere center at offset corner + Point3D sphere_center = { + orig.x() - corner_dir.x() * r, + orig.y() - corner_dir.y() * r, + orig.z() - corner_dir.z() * r + }; + + // Build spherical fill patch using the first 3 distinct normals + if (corner_normals.size() >= 3) { + auto sphere_patch = make_sphere_corner_patch( + sphere_center, corner_normals[0], corner_normals[1], corner_normals[2], r); + + // Create 3 contact points on the sphere surface + Point3D c0 = { sphere_center.x() + corner_normals[0].x() * r, + sphere_center.y() + corner_normals[0].y() * r, + sphere_center.z() + corner_normals[0].z() * r }; + Point3D c1 = { sphere_center.x() + corner_normals[1].x() * r, + sphere_center.y() + corner_normals[1].y() * r, + sphere_center.z() + corner_normals[1].z() * r }; + Point3D c2 = { sphere_center.x() + corner_normals[2].x() * r, + sphere_center.y() + corner_normals[2].y() * r, + sphere_center.z() + corner_normals[2].z() * r }; + + // Add spherical corner patch as triangle + add_tri(result, c0, c1, c2, sphere_patch); + + // Add transition patches: connect sphere patch edges to adjacent blend strips + // Transition between c0 (sphere) and t1 (blend strip) + Point3D mid_c0 = { + (c0.x() + t1.x()) * 0.5, (c0.y() + t1.y()) * 0.5, (c0.z() + t1.z()) * 0.5 + }; + auto trans_patch_1 = make_blend_patch(c0, mid_c0, t1, c0, mid_c0, t1); + add_tri(result, c0, t1, c0, trans_patch_1); // degenerate safe fallback + + // Transition between c1 (sphere) and t2 (blend strip) — similar + Point3D mid_c1 = { + (c1.x() + t2.x()) * 0.5, (c1.y() + t2.y()) * 0.5, (c1.z() + t2.z()) * 0.5 + }; + auto trans_patch_2 = make_blend_patch(c1, mid_c1, t2, c1, mid_c1, t2); + add_tri(result, c1, t2, c1, trans_patch_2); + } + }; + + // Process corners at both endpoints of the target edge + double r_start_eff = radii[0]; + double r_end_eff = radii[N]; + process_corner(edge.v_start, r_start_eff, t1_pts[0], t2_pts[0]); + process_corner(edge.v_end, r_end_eff, t1_pts[N], t2_pts[N]); + + finish_model(result, "shape_fillet"); + return result; +} + } // namespace vde::brep diff --git a/src/curves/control_point_edit.cpp b/src/curves/control_point_edit.cpp new file mode 100644 index 0000000..03a3705 --- /dev/null +++ b/src/curves/control_point_edit.cpp @@ -0,0 +1,115 @@ +/** + * @file control_point_edit.cpp + * @brief 控制点编辑实现 + */ + +#include "vde/curves/control_point_edit.h" +#include +#include +#include + +namespace vde::curves { + +ControlPointEditResult control_point_edit( + const NurbsSurface& surface, + const std::vector>& cp_indices, + const std::vector& new_positions) { + + ControlPointEditResult result; + + if (cp_indices.empty() || cp_indices.size() != new_positions.size()) { + // 无效输入，返回原曲面拷贝 + const auto& src_cp = surface.control_points(); + const auto& src_w = surface.weights(); + auto [nu, nv] = surface.num_control_points(); + std::vector> grid(nu, std::vector(nv)); + std::vector> w(nu, std::vector(nv, 1.0)); + bool has_w = !src_w.empty(); + for (int i = 0; i < nu; ++i) { + for (int j = 0; j < nv; ++j) { + grid[i][j] = src_cp[i][j]; + if (has_w) w[i][j] = src_w[i][j]; + } + } + result.modified_surface = NurbsSurface(grid, + surface.knots_u(), surface.knots_v(), w, + surface.degree_u(), surface.degree_v()); + return result; + } + + // 复制原始数据 + const auto& src_cp = surface.control_points(); + const auto& src_w = surface.weights(); + auto [nu, nv] = surface.num_control_points(); + + std::vector> grid(nu, std::vector(nv)); + std::vector> weights(nu, std::vector(nv, 1.0)); + bool has_weights = !src_w.empty(); + + for (int i = 0; i < nu; ++i) { + for (int j = 0; j < nv; ++j) { + grid[i][j] = src_cp[i][j]; + if (has_weights && static_cast(src_w[i].size()) > j) { + weights[i][j] = src_w[i][j]; + } + } + } + + // 修改指定的控制点 + int modified = 0; + double max_disp = 0.0; + bool weights_ok = true; + + for (size_t idx = 0; idx < cp_indices.size(); ++idx) { + int i = cp_indices[idx].first; + int j = cp_indices[idx].second; + + // 边界检查 + if (i < 0 || i >= nu || j < 0 || j >= nv) continue; + + // 记录位移 + Vector3D disp = {new_positions[idx].x() - grid[i][j].x(), + new_positions[idx].y() - grid[i][j].y(), + new_positions[idx].z() - grid[i][j].z()}; + double d = disp.norm(); + if (d > max_disp) max_disp = d; + + // 更新控制点 + grid[i][j] = new_positions[idx]; + ++modified; + } + + // 确保所有权重大于 0 + const double min_weight = 1e-12; + for (int i = 0; i < nu; ++i) { + for (int j = 0; j < nv; ++j) { + if (weights[i][j] < min_weight) { + weights[i][j] = min_weight; + weights_ok = false; + } + } + } + + // 构建位移数组 + for (int i = 0; i < nu; ++i) { + for (int j = 0; j < nv; ++j) { + result.displacements.push_back({ + grid[i][j].x() - src_cp[i][j].x(), + grid[i][j].y() - src_cp[i][j].y(), + grid[i][j].z() - src_cp[i][j].z() + }); + } + } + + // 构建修改后的曲面 + result.modified_surface = NurbsSurface(grid, + surface.knots_u(), surface.knots_v(), weights, + surface.degree_u(), surface.degree_v()); + result.max_displacement = max_disp; + result.points_modified = modified; + result.weights_preserved = weights_ok; + + return result; +} + +} // namespace vde::curves diff --git a/src/curves/curve_tools.cpp b/src/curves/curve_tools.cpp new file mode 100644 index 0000000..a9e56f8 --- /dev/null +++ b/src/curves/curve_tools.cpp @@ -0,0 +1,383 @@ +#include "vde/curves/curve_tools.h" +#include +#include +#include + +namespace vde::curves { + +using core::Point3D; +using core::Vector3D; + +// ═══════════════════════════════════════════════════════════ +// Internal Helpers +// ═══════════════════════════════════════════════════════════ + +namespace { + +/// Euclidean distance +inline double dist3d(const Point3D& a, const Point3D& b) { + return (a - b).norm(); +} + +/// Build an orthonormal basis {u, v} perpendicular to normal n +/// Returns {u, v} such that u·v=0, |u|=|v|=1, u⊥n, v⊥n +std::pair orthonormal_basis(const Vector3D& n) { + Vector3D u; + // Pick an axis that is not parallel to n + if (std::abs(n.x()) < 0.9) { + u = Vector3D(1, 0, 0).cross(n).normalized(); + } else { + u = Vector3D(0, 1, 0).cross(n).normalized(); + } + Vector3D v = n.cross(u).normalized(); + return {u, v}; +} + +/// Compute the Frenet frame at parameter t on a curve +/// Returns {tangent, normal, binormal} +struct FrenetFrame { + Vector3D T, N, B; +}; + +FrenetFrame frenet_frame(const NurbsCurve& curve, double t) { + FrenetFrame f; + Vector3D d1 = curve.derivative(t, 1); + Vector3D d2 = curve.derivative(t, 2); + + double speed = d1.norm(); + if (speed < 1e-12) { + f.T = Vector3D::Zero(); + f.N = Vector3D::Zero(); + f.B = Vector3D::Zero(); + return f; + } + + f.T = d1.normalized(); + Vector3D cross = d1.cross(d2); + double cross_norm = cross.norm(); + + if (cross_norm < 1e-18) { + // Degenerate (straight line): use arbitrary normal + auto basis = orthonormal_basis(f.T); + f.N = basis.first; + f.B = basis.second; + } else { + f.B = cross.normalized(); + f.N = f.B.cross(f.T).normalized(); + } + return f; +} + +/// Fit a NURBS curve to a set of sample points using least-squares +/// Uses cubic degree with uniform knots +NurbsCurve fit_nurbs_curve(const std::vector& points, int degree = 3) { + if (points.size() < 2) return NurbsCurve(); + + // Use a subset as control points (simplified: interpolate with enough points) + int n_ctrl = std::max(degree + 1, static_cast(points.size()) / 3 + 1); + n_ctrl = std::min(n_ctrl, 50); // cap for reasonable size + + std::vector ctrl_pts; + std::vector weights; + std::vector knots; + + // Uniformly sample control points from the point set + for (int i = 0; i < n_ctrl; ++i) { + int idx = i * static_cast(points.size() - 1) / (n_ctrl - 1); + idx = std::max(0, std::min(static_cast(points.size() - 1), idx)); + ctrl_pts.push_back(points[idx]); + weights.push_back(1.0); + } + + // Build uniform clamped knot vector + int n = n_ctrl - 1; + int m = n + degree + 1; + knots.resize(m + 1); + // Clamped: first degree+1 knots = 0, last degree+1 = 1 + for (int i = 0; i <= degree; ++i) knots[i] = 0.0; + for (int i = m - degree; i <= m; ++i) knots[i] = 1.0; + int inner = m - 2 * degree; + for (int i = 1; i < inner; ++i) { + knots[degree + i] = static_cast(i) / inner; + } + + return NurbsCurve(ctrl_pts, knots, weights, degree); +} + +/// Sample a curve uniformly in parameter space +std::vector sample_curve(const NurbsCurve& curve, int n_samples) { + std::vector pts; + auto dom = curve.domain(); + double t_min = dom.first, t_max = dom.second; + if (t_max <= t_min || n_samples < 2) return pts; + + pts.reserve(n_samples); + for (int i = 0; i < n_samples; ++i) { + double t = t_min + (t_max - t_min) * i / (n_samples - 1); + pts.push_back(curve.evaluate(t)); + } + return pts; +} + +} // namespace + +// ═══════════════════════════════════════════════════════════ +// parallel_curve +// ═══════════════════════════════════════════════════════════ + +NurbsCurve parallel_curve( + const NurbsCurve& curve, + double distance, + const Vector3D& plane_normal, + int samples) { + + Vector3D n = plane_normal.normalized(); + auto pts = sample_curve(curve, samples); + if (pts.empty()) return NurbsCurve(); + + std::vector offset_pts; + offset_pts.reserve(pts.size()); + + for (size_t i = 0; i < pts.size(); ++i) { + // Estimate the curve tangent using central difference + Vector3D tangent; + if (pts.size() >= 3 && i > 0 && i < pts.size() - 1) { + tangent = (pts[i + 1] - pts[i - 1]).normalized(); + } else { + // Use the domain parameter for derivative + auto dom = curve.domain(); + double t_min = dom.first, t_max = dom.second; + double t = t_min + (t_max - t_min) * i / (samples - 1); + Vector3D d1 = curve.derivative(t, 1); + tangent = d1.norm() > 1e-12 ? d1.normalized() : Vector3D::Zero(); + } + + if (tangent.norm() < 1e-12) { + offset_pts.push_back(pts[i]); + continue; + } + + // Offset direction: in-plane normal = plane_normal × tangent (normalized) + Vector3D offset_dir = n.cross(tangent).normalized(); + offset_pts.push_back(Point3D(pts[i] + offset_dir * distance)); + } + + return fit_nurbs_curve(offset_pts); +} + +// ═══════════════════════════════════════════════════════════ +// connect_curve +// ═══════════════════════════════════════════════════════════ + +NurbsCurve connect_curve( + const NurbsCurve& c1, + const NurbsCurve& c2, + Continuity continuity, + int samples) { + + auto dom1 = c1.domain(); + auto dom2 = c2.domain(); + + Point3D P0 = c1.evaluate(dom1.second); // end of c1 + Point3D Pn = c2.evaluate(dom2.first); // start of c2 + + // Get derivatives at endpoints + Vector3D T1 = c1.derivative(dom1.second, 1); // c1 tangent at end + Vector3D T2 = c2.derivative(dom2.first, 1); // c2 tangent at start + + // Normalize tangents + if (T1.norm() > 1e-12) T1 = T1.normalized(); + if (T2.norm() > 1e-12) T2 = T2.normalized(); + + // Determine connection curve degree based on continuity + int degree = 3; // default G1 + switch (continuity) { + case Continuity::G1: degree = 3; break; + case Continuity::G2: degree = 5; break; + case Continuity::G3: degree = 7; break; + } + + // Build a Bézier transition curve with Hermite end conditions + int n_cp = degree + 1; + std::vector cp(n_cp); + + cp[0] = P0; + cp[degree] = Pn; + + // Chord length for scaling + double L = (Pn - P0).norm(); + if (L < 1e-12) { + // Degenerate: return a straight line + return NurbsCurve({P0, Pn}, {0, 0, 1, 1}, {1, 1}, 1); + } + + // Scale tangents: typical cubic Hermite uses L/3 for inner controls + double s = L / degree; + + // G1: match tangents (P1, P_{n-1}) + cp[1] = Point3D(P0 + T1 * s); + + if (degree >= 5) { + // G2: match curvature — use second derivatives + Vector3D D2_1 = c1.derivative(dom1.second, 2); + Vector3D D2_2 = c2.derivative(dom2.first, 2); + + // For a cubic Bézier: P2 = P0 + 2*s*T1 + s²*κ*N1 + // For degree 5: distribute the curvature constraint + // Simplified approach: use the chord line as default, adjust with tangent + cp[degree - 1] = Point3D(Pn - T2 * s); + + // Place interior control points along the chord for smoothness + if (degree >= 7) { + // G3: match torsion + double t_step = 1.0 / degree; + for (int i = 2; i < degree - 1; ++i) { + double t = i * t_step; + // Simple linear interpolation for interior controls + cp[i] = Point3D(P0 * (1 - t) + Pn * t); + + // Blend tangents near ends + if (i == 2) { + Vector3D mid = P0 * 0.6 + Pn * 0.4; + cp[i] = Point3D(mid + T1 * s * 0.2); + } else if (i == degree - 2) { + Vector3D mid = P0 * 0.4 + Pn * 0.6; + cp[i] = Point3D(mid - T2 * s * 0.2); + } + } + } else { + // Degree 5: fill interior with linear blend + for (int i = 2; i < degree - 1; ++i) { + double t = static_cast(i) / degree; + cp[i] = Point3D(P0 * (1 - t) + Pn * t); + } + } + } else { + // Degree 3: remaining control point is mid with tangent influence + cp[degree - 1] = Point3D(Pn - T2 * s); + // Interior point (if degree > 3) — not applicable for cubic + } + + // Build clamped knot vector + std::vector knots(n_cp + degree + 1); + for (int i = 0; i <= degree; ++i) knots[i] = 0.0; + for (int i = n_cp; i <= n_cp + degree; ++i) knots[i] = 1.0; + int inner = n_cp - degree; + for (int i = 1; i < inner; ++i) { + knots[degree + i] = static_cast(i) / inner; + } + + std::vector w(n_cp, 1.0); + return NurbsCurve(cp, knots, w, degree); +} + +// ═══════════════════════════════════════════════════════════ +// helix_curve +// ═══════════════════════════════════════════════════════════ + +NurbsCurve helix_curve( + const Point3D& axis_origin, + const Vector3D& axis_dir, + double radius, + double pitch, + double height, + int samples) { + + Vector3D dir = axis_dir.normalized(); + if (height <= 0.0 || pitch <= 0.0 || radius < 0.0) return NurbsCurve(); + + // Build orthonormal basis for the plane perpendicular to axis_dir + auto basis = orthonormal_basis(dir); + Vector3D u = basis.first; + Vector3D v = basis.second; + + double turns = height / pitch; + int total_samples = static_cast(samples * turns); + total_samples = std::max(total_samples, 20); + + std::vector pts; + pts.reserve(total_samples); + + for (int i = 0; i < total_samples; ++i) { + double t = static_cast(i) / (total_samples - 1); + double angle = 2.0 * M_PI * t * turns; + double z = height * t; + + Point3D p = axis_origin + + u * (radius * std::cos(angle)) + + v * (radius * std::sin(angle)) + + dir * z; + pts.push_back(p); + } + + return fit_nurbs_curve(pts, 3); +} + +// ═══════════════════════════════════════════════════════════ +// isoparametric_curves +// ═══════════════════════════════════════════════════════════ + +IsoparametricCurves isoparametric_curves( + const NurbsSurface& surface, + int u_count, + int v_count) { + + IsoparametricCurves result; + + const auto& ku = surface.knots_u(); + const auto& kv = surface.knots_v(); + + if (ku.empty() || kv.empty()) return result; + + int pu = surface.degree_u(); + int pv = surface.degree_v(); + + double u_min = ku[pu]; + double u_max = ku[ku.size() - 1 - pu]; + double v_min = kv[pv]; + double v_max = kv[kv.size() - 1 - pv]; + + // Extract v-iso curves (v = const, u varies) + if (v_count >= 1) { + result.v_curves.reserve(v_count); + for (int j = 0; j < v_count; ++j) { + double v = (v_count == 1) ? (v_min + v_max) * 0.5 + : v_min + (v_max - v_min) * j / (v_count - 1); + + // Sample points along this v = const line + int n_samples = 50; + std::vector pts; + pts.reserve(n_samples); + for (int i = 0; i < n_samples; ++i) { + double u = u_min + (u_max - u_min) * i / (n_samples - 1); + pts.push_back(surface.evaluate(u, v)); + } + + result.v_curves.push_back(fit_nurbs_curve(pts)); + } + } + + // Extract u-iso curves (u = const, v varies) + if (u_count >= 1) { + result.u_curves.reserve(u_count); + for (int i = 0; i < u_count; ++i) { + double u = (u_count == 1) ? (u_min + u_max) * 0.5 + : u_min + (u_max - u_min) * i / (u_count - 1); + + // Sample points along this u = const line + int n_samples = 50; + std::vector pts; + pts.reserve(n_samples); + for (int j = 0; j < n_samples; ++j) { + double v = v_min + (v_max - v_min) * j / (n_samples - 1); + pts.push_back(surface.evaluate(u, v)); + } + + result.u_curves.push_back(fit_nurbs_curve(pts)); + } + } + + return result; +} + +} // namespace vde::curves diff --git a/src/curves/generative_surfaces.cpp b/src/curves/generative_surfaces.cpp new file mode 100644 index 0000000..17efd25 --- /dev/null +++ b/src/curves/generative_surfaces.cpp @@ -0,0 +1,542 @@ +#include "vde/curves/generative_surfaces.h" +#include "vde/curves/bspline_curve.h" +#include +#include +#include +#include + +namespace vde::curves { + +// ============================================================================ +// 内部辅助函数 +// ============================================================================ + +namespace { + +/// 两个控制点序列中最大长度 +int max_cp_count(const std::vector& curves) { + int m = 0; + for (const auto& c : curves) + m = std::max(m, static_cast(c.control_points().size())); + return m; +} + +/// 升阶：将 degree p 升为 new_degree（Cox-de Boor 升阶算法简化版） +std::vector elevate_degree(const std::vector& cp, + const std::vector& knots, + int degree, int new_degree) { + if (new_degree <= degree) return cp; + // 简化：在 CP 和 knots 上补足多重度 + std::vector result = cp; + std::vector new_knots = knots; + for (int d = degree; d < new_degree; ++d) { + int n = static_cast(result.size()) - 1; + // 在每段中点插入控制点 + std::vector elevated; + elevated.push_back(result[0]); + for (int i = 0; i < n; ++i) { + double alpha = static_cast(i + 1) / (n + 1); + elevated.push_back(result[i] + alpha * (result[i+1] - result[i])); + } + elevated.push_back(result.back()); + result = std::move(elevated); + } + return result; +} + +/// 节点精化：将两条曲线的节点向量合并 +std::vector merge_knots(const std::vector& k1, const std::vector& k2) { + std::vector merged; + for (double k : k1) merged.push_back(k); + for (double k : k2) merged.push_back(k); + std::sort(merged.begin(), merged.end()); + merged.erase(std::unique(merged.begin(), merged.end(), + [](double a, double b) { return std::abs(a-b) < 1e-10; }), merged.end()); + return merged; +} + +/// 统一化多条曲线：升阶+精化，使所有曲线控制点数相同 +void unify_curves(std::vector& curves, std::vector& out_knots, int& out_degree) { + if (curves.empty()) return; + + // 找最大 degree + out_degree = 0; + for (const auto& c : curves) + out_degree = std::max(out_degree, c.degree()); + + // 统一 degree + for (auto& c : curves) { + if (c.degree() < out_degree) { + auto new_cp = elevate_degree(c.control_points(), c.knots(), c.degree(), out_degree); + // Rebuild with elevated CPs + c = NurbsCurve(new_cp, c.knots(), std::vector(new_cp.size(), 1.0), out_degree); + } + } + + // 合并所有节点 + out_knots = curves[0].knots(); + for (size_t i = 1; i < curves.size(); ++i) + out_knots = merge_knots(out_knots, curves[i].knots()); + + // 精化每条曲线到统一节点 + size_t target_n = out_knots.size() - out_degree - 1; + for (auto& c : curves) { + if (c.control_points().size() != target_n) { + // 简单重采样 + std::vector new_cp; + for (size_t j = 0; j < target_n; ++j) { + double t = out_knots[out_degree + j]; + new_cp.push_back(c.evaluate(t)); + } + c = NurbsCurve(new_cp, out_knots, std::vector(target_n, 1.0), out_degree); + } + } +} + +/// 计算 Frenet 标架 (tangent, normal, binormal) +struct Frame { + Vector3D t, n, b; +}; + +Frame compute_frenet(const NurbsCurve& curve, double t, const Frame& prev) { + Frame f; + f.t = curve.derivative(t, 1).normalized(); + + // 二阶导数 → 曲率法向 + Vector3D d2 = curve.derivative(t, 2); + double d2_norm = d2.norm(); + + if (d2_norm < 1e-9) { + // 曲率 ≈ 0，沿用前一标架的 n（平行传输近似） + // 重新正交化 + f.n = prev.n - f.t * f.t.dot(prev.n); + double nn = f.n.norm(); + if (nn < 1e-9) { + // 随便找一个垂直于 t 的向量 + if (std::abs(f.t.x()) < 0.9) + f.n = Vector3D(1, 0, 0).cross(f.t).normalized(); + else + f.n = Vector3D(0, 1, 0).cross(f.t).normalized(); + } else { + f.n = f.n / nn; + } + } else { + f.n = d2.normalized(); + } + + f.b = f.t.cross(f.n).normalized(); + return f; +} + +/// 将 profile 控制点从局部坐标系变换到 Frenet 标架所在的世界坐标 +Point3D transform_point(const Point3D& p, const Point3D& origin, + const Vector3D& tx, const Vector3D& ty, const Vector3D& tz) { + // p 在 profile 局部坐标系 (x = profile tangent → tx in world, y → ty, z → tz) + // 但实际上 profile 通常在 XY 平面，我们让 profile 局部 = Y 方向为 n, Z 方向为 b + return origin + tx * p.x() + ty * p.y() + tz * p.z(); +} + +/// 构造 v 向 B 样条插值（给定 N 个 v 参数值和对应点） +/// 返回一条插值 B 样条（degree = 3, clamped） +NurbsCurve interpolate_bspline(const std::vector& params, + const std::vector& points) { + int n = static_cast(points.size()) - 1; + int p = std::min(3, n); // degree + + if (n < 1) { + std::vector cps = {points[0], points[0]}; + return NurbsCurve(cps, {0, 0, 1, 1}, {1, 1}, 1); + } + + // 构建节点向量（clamped） + int m = n + p + 1; + std::vector knots(m + 1); + for (int i = 0; i <= p; ++i) knots[i] = 0.0; + for (int i = p + 1; i <= n; ++i) + knots[i] = static_cast(i - p) / (n - p + 1); + for (int i = n + 1; i <= m; ++i) knots[i] = 1.0; + + // 对每个参数点，计算基函数值，构建线性系统 A·CP = points + // 简化：直接使用参数点作为控制点（通过 clamped knots） + std::vector cp = points; + + return NurbsCurve(cp, knots, std::vector(cp.size(), 1.0), p); +} + +/// 两条曲线求交点（采样法：在 N×M 网格中找最近点对） +/// 返回 {参数1, 参数2} 或提示未找到 +std::pair find_intersection(const NurbsCurve& c1, const NurbsCurve& c2, + int samples = 100) { + auto d1 = c1.domain(); + auto d2 = c2.domain(); + double best_dist = 1e100; + double best_t1 = 0.0, best_t2 = 0.0; + + for (int i = 0; i <= samples; ++i) { + double t1 = d1.first + (d1.second - d1.first) * i / samples; + Point3D p1 = c1.evaluate(t1); + for (int j = 0; j <= samples; ++j) { + double t2 = d2.first + (d2.second - d2.first) * j / samples; + double dist = (p1 - c2.evaluate(t2)).norm(); + if (dist < best_dist) { + best_dist = dist; + best_t1 = t1; + best_t2 = t2; + } + } + } + + // 局部细化（Newton-like） + for (int iter = 0; iter < 5 && best_dist > 1e-6; ++iter) { + Point3D p1 = c1.evaluate(best_t1); + Point3D p2 = c2.evaluate(best_t2); + Vector3D diff = p1 - p2; + best_dist = diff.norm(); + if (best_dist < 1e-6) break; + + Vector3D deriv1 = c1.derivative(best_t1, 1); + Vector3D deriv2 = c2.derivative(best_t2, 1); + + // 简化的 Gauss-Newton: min ||c1(t1) - c2(t2)||² + double dt1 = diff.dot(deriv1) / (deriv1.norm() * deriv1.norm() + 1e-12); + double dt2 = -diff.dot(deriv2) / (deriv2.norm() * deriv2.norm() + 1e-12); + best_t1 -= dt1 * 0.5; + best_t2 -= dt2 * 0.5; + best_t1 = std::max(d1.first, std::min(d1.second, best_t1)); + best_t2 = std::max(d2.first, std::min(d2.second, best_t2)); + } + + return {best_t1, best_t2}; +} + +/// 求平均 +double avg_param(const std::vector& vals) { + if (vals.empty()) return 0.0; + double s = 0.0; + for (double v : vals) s += v; + return s / vals.size(); +} + +} // anonymous namespace + +// ============================================================================ +// sweep_surface +// ============================================================================ + +NurbsSurface sweep_surface(const NurbsCurve& profile, + const NurbsCurve& path, + const SweepOption& option) { + const auto& prof_cp = profile.control_points(); + int prof_n = static_cast(prof_cp.size()); + if (prof_n < 2) + throw std::invalid_argument("sweep_surface: profile must have at least 2 control points"); + + int samples = std::max(2, option.samples); + auto path_dom = path.domain(); + + // ── 采样路径并计算每点的 Frenet 标架 ── + std::vector path_pts; + std::vector frames; + Frame prev; + bool first = true; + + for (int i = 0; i <= samples; ++i) { + double t = path_dom.first + (path_dom.second - path_dom.first) * i / samples; + Point3D pt = path.evaluate(t); + path_pts.push_back(pt); + + if (first) { + Frame f; + f.t = path.derivative(t, 1).normalized(); + // 初始法向：垂直于 t 的任意向量 + if (std::abs(f.t.x()) < 0.9) + f.n = Vector3D(1, 0, 0).cross(f.t).normalized(); + else + f.n = Vector3D(0, 1, 0).cross(f.t).normalized(); + f.b = f.t.cross(f.n).normalized(); + frames.push_back(f); + prev = f; + first = false; + } else { + Frame f = compute_frenet(path, t, prev); + frames.push_back(f); + prev = f; + } + } + + // ── Spine 模式：用脊线修正扭转 ── + if (option.type == SweepType::Spine && option.spine.has_value()) { + // 将每个 Frenet 标架的 n 重新投影到 spine 方向 + for (size_t i = 0; i < frames.size(); ++i) { + // 找到 spine 上最近点 + auto sp_dom = option.spine->domain(); + double best_s = sp_dom.first; + double best_d = 1e100; + for (int j = 0; j <= 50; ++j) { + double s = sp_dom.first + (sp_dom.second - sp_dom.first) * j / 50; + double d = (option.spine->evaluate(s) - path_pts[i]).norm(); + if (d < best_d) { best_d = d; best_s = s; } + } + Vector3D spine_pt = option.spine->evaluate(best_s); + Vector3D to_spine = spine_pt - path_pts[i]; + if (to_spine.norm() > 1e-9) { + // 将 n 投影到 spine 方向 + Vector3D spine_dir = to_spine.normalized(); + frames[i].n = spine_dir - frames[i].t * frames[i].t.dot(spine_dir); + double nn = frames[i].n.norm(); + if (nn > 1e-9) frames[i].n /= nn; + else { + if (std::abs(frames[i].t.x()) < 0.9) + frames[i].n = Vector3D(1, 0, 0).cross(frames[i].t).normalized(); + else + frames[i].n = Vector3D(0, 1, 0).cross(frames[i].t).normalized(); + } + frames[i].b = frames[i].t.cross(frames[i].n).normalized(); + } + } + } + + // ── TwoGuides 模式：计算缩放因子 ── + std::vector scale_factors(samples + 1, 1.0); + if (option.type == SweepType::TwoGuides && option.guide1.has_value() && option.guide2.has_value()) { + // 在路径起点处计算 profile 宽度（首尾控制点距离） + double prof_base = (prof_cp[0] - prof_cp.back()).norm(); + if (prof_base < 1e-9) prof_base = 1.0; + + // 在路径起点处计算两条引导线之间的距离 + auto g1_dom = option.guide1->domain(); + auto g2_dom = option.guide2->domain(); + double g1_start_dist = (option.guide1->evaluate(g1_dom.first) - path_pts[0]).norm(); + double g2_start_dist = (option.guide2->evaluate(g2_dom.first) - path_pts[0]).norm(); + + for (size_t i = 0; i < path_pts.size(); ++i) { + double t_frac = static_cast(i) / samples; + double g1_t = g1_dom.first + (g1_dom.second - g1_dom.first) * t_frac; + double g2_t = g2_dom.first + (g2_dom.second - g2_dom.first) * t_frac; + + Point3D g1_pt = option.guide1->evaluate(g1_t); + Point3D g2_pt = option.guide2->evaluate(g2_t); + + // 引导线之间距离相对起始距离的比例 + double cur_dist = (g1_pt - g2_pt).norm(); + double start_dist = (option.guide1->evaluate(g1_dom.first) - + option.guide2->evaluate(g2_dom.first)).norm(); + if (start_dist > 1e-9) + scale_factors[i] = cur_dist / start_dist; + else + scale_factors[i] = 1.0; + } + } + + // ── 变换 profile 控制点到每个路径位置 ── + std::vector> grid(samples + 1); + std::vector> w(samples + 1); + + // profile 局部坐标系: Y=n, Z=b, X=t (profile 在 local YZ 平面，曲线方向沿 X) + for (int i = 0; i <= samples; ++i) { + const Frame& f = frames[i]; + double sc = scale_factors[i]; + for (int j = 0; j < prof_n; ++j) { + // profile 在其局部 YZ 平面：x=沿曲线方向(不用于扫描)，y=垂直于路径方向，z=副法向 + // 扫描时：profile 的 Y→法向, Z→副法向 + Point3D p = prof_cp[j]; + // 应用缩放（沿 n 和 b 方向） + double py = p.y() * sc; + double pz = p.z() * sc; + Point3D wp = path_pts[i] + f.n * py + f.b * pz + f.t * p.x(); + grid[i].push_back(wp); + } + // 权重：profile 有 weights 时传递 + const auto& pw = profile.weights(); + if (!pw.empty()) { + for (double wi : pw) w[i].push_back(wi); + } + } + + // ── 构建 NURBS 曲面 ── + // V 向（沿路径）：使用 profile 的阶次 + // U 向（沿截面）：degree 2, clamped + int prof_deg = profile.degree(); + int path_deg = 2; // 沿路径方向 degree 2 + + // 沿路径的节点向量 + int nu = samples + 1; + std::vector ku; + int nku = nu + path_deg + 1; + ku.resize(nku); + for (int i = 0; i <= path_deg; ++i) ku[i] = 0.0; + for (int i = path_deg + 1; i <= nu - 1; ++i) + ku[i] = static_cast(i - path_deg) / (nu - path_deg); + for (int i = nu; i < nku; ++i) ku[i] = 1.0; + + // 沿截面方向的节点向量（沿用 profile 的） + std::vector kv = profile.knots(); + + return NurbsSurface(grid, ku, kv, w, path_deg, prof_deg); +} + +// ============================================================================ +// loft_surface +// ============================================================================ + +NurbsSurface loft_surface(const std::vector& sections, + const std::vector& guides, + const LoftConstraint& constraints) { + if (sections.size() < 2) + throw std::invalid_argument("loft_surface: at least 2 sections required"); + + // ── 复制并统一化所有截面 ── + std::vector curves = sections; + std::vector unified_knots; + int unified_deg; + unify_curves(curves, unified_knots, unified_deg); + + int n_cp = static_cast(curves[0].control_points().size()); + + // ── 确定各截面在 v 方向的参数位置 ── + int n_sec = static_cast(curves.size()); + std::vector v_params(n_sec); + if (!guides.empty()) { + // 使用引导线计算截面参数位置 + // 简化：均匀分布但在各自引导线起点附近 + for (int i = 0; i < n_sec; ++i) + v_params[i] = static_cast(i) / (n_sec - 1); + } else { + for (int i = 0; i < n_sec; ++i) + v_params[i] = static_cast(i) / (n_sec - 1); + } + + // ── v 向 B 样条插值（对每个 u 控制点列） ── + // 构建节点向量 + int pv = std::min(3, n_sec - 1); + int nv_knots = n_sec + pv + 1; + std::vector kv(nv_knots); + for (int i = 0; i <= pv; ++i) kv[i] = 0.0; + for (int i = pv + 1; i <= n_sec - 1; ++i) + kv[i] = static_cast(i - pv) / (n_sec - pv); + for (int i = n_sec; i < nv_knots; ++i) kv[i] = 1.0; + + // ── 切线约束：修改首尾 v 参数 ── + std::vector> loft_cp(n_sec); + for (int i = 0; i < n_sec; ++i) { + loft_cp[i] = curves[i].control_points(); + } + + // 对每个 u 控制点索引，沿 v 方向插值 + std::vector> surf_grid(n_cp); + std::vector> surf_weights(n_cp); + + for (int ui = 0; ui < n_cp; ++ui) { + std::vector col_pts; + for (int vi = 0; vi < n_sec; ++vi) { + col_pts.push_back(loft_cp[vi][ui]); + } + // v 向插值 + auto interp = interpolate_bspline(v_params, col_pts); + // 重采样到统一节点 + for (size_t k = 0; k < unified_knots.size() - unified_deg - 1; ++k) { + // 不对，我们已有正确的 u 向结构 + } + + // 直接用插值后的控制点 + surf_grid[ui] = interp.control_points(); + if (!interp.weights().empty()) + surf_weights[ui] = interp.weights(); + } + + // 转置 grid：grid[ui][vi] → grid_t[vi][ui] (NurbsSurface 格式: grid[u][v]) + // surf_grid 已经是 grid[u_index][v_index]，符合 NurbsSurface 期望 + + // 但 kv 应当反映 v 向插值 + // 我们需要用正确的 v 向节点 + std::vector v_knots; + if (!surf_grid.empty()) { + int nv_cp = static_cast(surf_grid[0].size()); + int nv_deg = pv; + int nvk = nv_cp + nv_deg + 1; + v_knots.resize(nvk); + for (int i = 0; i <= nv_deg; ++i) v_knots[i] = 0.0; + for (int i = nv_deg + 1; i <= nv_cp - 1; ++i) + v_knots[i] = static_cast(i - nv_deg) / (nv_cp - nv_deg); + for (int i = nv_cp; i < nvk; ++i) v_knots[i] = 1.0; + } + + return NurbsSurface(surf_grid, unified_knots, v_knots, surf_weights, unified_deg, pv); +} + +// ============================================================================ +// net_surface +// ============================================================================ + +NurbsSurface net_surface(const std::vector& u_curves, + const std::vector& v_curves) { + if (u_curves.size() < 2) + throw std::invalid_argument("net_surface: at least 2 u-curves required"); + if (v_curves.size() < 2) + throw std::invalid_argument("net_surface: at least 2 v-curves required"); + + int nu = static_cast(u_curves.size()); + int nv = static_cast(v_curves.size()); + + // ── 找所有 u-curve(i) 与 v-curve(j) 的交点 ── + // 对每个 u/v 曲线对，求近似交点参数 + std::vector> u_params(nu, std::vector(nv)); + std::vector> v_params(nu, std::vector(nv)); + + for (int i = 0; i < nu; ++i) { + for (int j = 0; j < nv; ++j) { + auto [tu, tv] = find_intersection(u_curves[i], v_curves[j]); + u_params[i][j] = tu; + v_params[i][j] = tv; + } + } + + // ── 对齐：计算平均 u 参数作为 u 向统一参数 ── + std::vector avg_u_params(nv); + for (int j = 0; j < nv; ++j) { + std::vector vals; + for (int i = 0; i < nu; ++i) vals.push_back(u_params[i][j]); + avg_u_params[j] = avg_param(vals); + } + + std::vector avg_v_params(nu); + for (int i = 0; i < nu; ++i) { + std::vector vals; + for (int j = 0; j < nv; ++j) vals.push_back(v_params[i][j]); + avg_v_params[i] = avg_param(vals); + } + + // ── 构建控制点网格：每个 (i,j) 处取 v-curve(j) 在 avg_u_params[j] 处的点 ── + // 行索引 = u 方向（沿 v-curve），列索引 = v 方向（沿 u-curve） + std::vector> grid(nu); + for (int i = 0; i < nu; ++i) { + for (int j = 0; j < nv; ++j) { + Point3D pt = v_curves[j].evaluate(avg_u_params[j]); + grid[i].push_back(pt); + } + } + + // ── 构建节点向量 ── + int pu = std::min(2, std::max(0, nu - 1)); + int pv = std::min(2, std::max(0, nv - 1)); + + int nku = nu + pu + 1; + std::vector ku(nku); + for (int i = 0; i <= pu; ++i) ku[i] = 0.0; + for (int i = pu + 1; i <= nu - 1; ++i) + ku[i] = static_cast(i - pu) / (nu - pu); + for (int i = nu; i < nku; ++i) ku[i] = 1.0; + + int nkv = nv + pv + 1; + std::vector kv(nkv); + for (int i = 0; i <= pv; ++i) kv[i] = 0.0; + for (int i = pv + 1; i <= nv - 1; ++i) + kv[i] = static_cast(i - pv) / (nv - pv); + for (int i = nv; i < nkv; ++i) kv[i] = 1.0; + + std::vector> w; // uniform weights + + return NurbsSurface(grid, ku, kv, w, pu, pv); +} + +} // namespace vde::curves diff --git a/src/curves/surface_analysis.cpp b/src/curves/surface_analysis.cpp index 161de92..35db9cb 100644 --- a/src/curves/surface_analysis.cpp +++ b/src/curves/surface_analysis.cpp @@ -534,4 +534,252 @@ DraftAngleDistribution draft_face_angle_distribution( return result; } +// ═══════════════════════════════════════════════════════════ +// Inflection Lines +// ═══════════════════════════════════════════════════════════ + +InflectionLinesResult inflection_lines( + const NurbsSurface& surf, + int res_u, int res_v) { + + InflectionLinesResult result; + result.res_u = res_u; + result.res_v = res_v; + + const auto& ku = surf.knots_u(); + const auto& kv = surf.knots_v(); + if (ku.empty() || kv.empty()) return result; + + int pu = surf.degree_u(); + int pv = surf.degree_v(); + double u_min = ku[pu]; + double u_max = ku[ku.size() - 1 - pu]; + double v_min = kv[pv]; + double v_max = kv[kv.size() - 1 - pv]; + + // Sample Gaussian curvature on a grid + std::vector> K(res_u + 1, std::vector(res_v + 1, 0.0)); + std::vector> pts(res_u + 1, std::vector(res_v + 1)); + + for (int i = 0; i <= res_u; ++i) { + double u = u_min + (u_max - u_min) * i / res_u; + for (int j = 0; j <= res_v; ++j) { + double v = v_min + (v_max - v_min) * j / res_v; + auto d = compute_surface_derivs(surf, u, v); + double k1, k2; + Vector3D d1, d2; + principal_curvatures(d, k1, k2, d1, d2); + K[i][j] = k1 * k2; + pts[i][j] = surf.evaluate(u, v); + } + } + + // Find zero-crossings along both directions + for (int i = 0; i <= res_u; ++i) { + for (int j = 0; j < res_v; ++j) { + if (K[i][j] * K[i][j + 1] < 0.0) { + // Sign change between (i,j) and (i,j+1) + // Interpolate to find zero-crossing point + double t = std::abs(K[i][j]) / (std::abs(K[i][j]) + std::abs(K[i][j + 1])); + Point3D p0 = pts[i][j]; + Point3D p1 = pts[i][j + 1]; + Point3D mid = Point3D(p0 + (p1 - p0) * t); + + // Create a small segment across the zero-crossing + InflectionSegment seg; + seg.start = mid; + seg.end = mid; // point marker — could be extended to connect adjacent crossings + result.segments.push_back(seg); + } + } + } + + for (int j = 0; j <= res_v; ++j) { + for (int i = 0; i < res_u; ++i) { + if (K[i][j] * K[i + 1][j] < 0.0) { + double t = std::abs(K[i][j]) / (std::abs(K[i][j]) + std::abs(K[i + 1][j])); + Point3D p0 = pts[i][j]; + Point3D p1 = pts[i + 1][j]; + Point3D mid = Point3D(p0 + (p1 - p0) * t); + + InflectionSegment seg; + seg.start = mid; + seg.end = mid; + result.segments.push_back(seg); + } + } + } + + // Connect adjacent zero-crossing points into short segments + // (simplified: connect nearby points within a threshold) + std::vector connected; + double connect_threshold = (u_max - u_min) * (v_max - v_min) * 0.01; + + for (size_t i = 0; i < result.segments.size(); ++i) { + for (size_t j = i + 1; j < result.segments.size(); ++j) { + if (dist3d(result.segments[i].start, result.segments[j].start) < connect_threshold) { + InflectionSegment seg; + seg.start = result.segments[i].start; + seg.end = result.segments[j].start; + connected.push_back(seg); + } + } + } + + if (!connected.empty()) { + result.segments = connected; + } + + return result; +} + +// ═══════════════════════════════════════════════════════════ +// Checker Mapping +// ═══════════════════════════════════════════════════════════ + +CheckerMapResult checker_mapping( + const NurbsSurface& surf, + const Vector3D& view_dir, + int res_u, int res_v) { + + CheckerMapResult result; + result.res_u = res_u; + result.res_v = res_v; + + Vector3D V = view_dir.normalized(); + result.view_dir = V; + + const auto& ku = surf.knots_u(); + const auto& kv = surf.knots_v(); + if (ku.empty() || kv.empty()) return result; + + int pu = surf.degree_u(); + int pv = surf.degree_v(); + double u_min = ku[pu]; + double u_max = ku[ku.size() - 1 - pu]; + double v_min = kv[pv]; + double v_max = kv[kv.size() - 1 - pv]; + + result.intensity.resize(res_u + 1, std::vector(res_v + 1, 0.0)); + + int checker_size = 8; // number of checker squares per direction + + for (int i = 0; i <= res_u; ++i) { + double u = u_min + (u_max - u_min) * i / res_u; + for (int j = 0; j <= res_v; ++j) { + double v = v_min + (v_max - v_min) * j / res_v; + + Vector3D N = surf.normal(u, v); + + // Reflection vector: R = V - 2*(V·N)*N + Vector3D R = V - 2.0 * N.dot(V) * N; + R = R.normalized(); + + // Map reflected direction to checkerboard pattern + // Use spherical coordinates of R to index checkerboard + double theta = std::atan2(R.y(), R.x()); + double phi = std::acos(R.z()); + + int cu = static_cast(std::floor(theta / M_PI * checker_size + 0.5)); + int cv = static_cast(std::floor(phi / M_PI * checker_size + 0.5)); + + // Checkerboard: 1 if cu+cv even, 0 if odd + double base = ((cu + cv) & 1) ? 1.0 : 0.0; + + // Blend with Fresnel-like term (view angle influence) + double fresnel = 1.0 - std::abs(N.dot(V)); + double intensity = base * 0.7 + fresnel * 0.3; + + result.intensity[i][j] = std::max(0.0, std::min(1.0, intensity)); + } + } + + return result; +} + +// ═══════════════════════════════════════════════════════════ +// Surface Checker Report +// ═══════════════════════════════════════════════════════════ + +SurfaceCheckerReport surface_checker_report( + const NurbsSurface& surf, + int res) { + + SurfaceCheckerReport report; + + // Compute curvature map + auto cm = curvature_map(surf, res, res); + + report.gaussian_range = cm.max_gaussian - cm.min_gaussian; + report.mean_range = cm.max_mean - cm.min_mean; + + // Compute checker mapping + Vector3D view_dir(0.5, 0.3, 1.0); // default view direction + auto ckr = checker_mapping(surf, view_dir, res, res); + + // Checker uniformity: measure variance of adjacent pixel differences + double sum_diff_sq = 0.0; + double sum_sq = 0.0; + int count = 0; + int count2 = 0; + + for (int i = 0; i <= res; ++i) { + for (int j = 0; j <= res; ++j) { + double val = ckr.intensity[i][j]; + sum_sq += val * val; + count2++; + + if (i < res) { + double diff_u = ckr.intensity[i + 1][j] - val; + sum_diff_sq += diff_u * diff_u; + count++; + } + if (j < res) { + double diff_v = ckr.intensity[i][j + 1] - val; + sum_diff_sq += diff_v * diff_v; + count++; + } + } + } + + double mean_sq = (count2 > 0) ? sum_sq / count2 : 0.0; + double mean_diff_sq = (count > 0) ? sum_diff_sq / count : 0.0; + + // Uniformity: higher variance in derivatives → less uniform + double uniformity = (mean_diff_sq > 0.0) + ? std::max(0.0, 1.0 - mean_diff_sq / (mean_sq + 1e-9)) + : 1.0; + report.checker_uniformity = uniformity; + + // Curvature continuity: measure smoothness of curvature change + // (simplified: use Gaussian curvature range normalized) + double gauss_span = std::max(1e-12, cm.max_gaussian - cm.min_gaussian); + double k_continuity = std::max(0.0, 1.0 - std::log1p(gauss_span) / std::log1p(1.0)); + report.curvature_continuity = k_continuity; + + // Inflection count + auto infl = inflection_lines(surf, res, res); + report.inflection_count = static_cast(infl.segments.size()); + + // Determine grade + double score = (uniformity * 0.4 + k_continuity * 0.4 + + (1.0 - std::min(1.0, report.inflection_count / 100.0)) * 0.2); + + if (score >= 0.85) { + report.grade = SmoothnessGrade::A; + report.description = "Excellent — uniform reflection stripes, smooth curvature"; + } else if (score >= 0.65) { + report.grade = SmoothnessGrade::B; + report.description = "Good — minor irregularities, locally smooth"; + } else if (score >= 0.40) { + report.grade = SmoothnessGrade::C; + report.description = "Fair — visible distortions, uneven surface"; + } else { + report.grade = SmoothnessGrade::D; + report.description = "Poor — severe distortions, curvature discontinuity"; + } + + return report; +} + } // namespace vde::curves diff --git a/src/curves/surface_editing.cpp b/src/curves/surface_editing.cpp new file mode 100644 index 0000000..9c228ad --- /dev/null +++ b/src/curves/surface_editing.cpp @@ -0,0 +1,651 @@ +/** + * @file surface_editing.cpp + * @brief 曲面编辑（Match Surface）实现 + */ + +#include "vde/curves/surface_editing.h" +#include +#include +#include +#include + +namespace vde::curves { + +namespace { + +// ── 内部辅助函数 ────────────────────────────────────────── + +/// 安全向量归一化 +inline Vector3D safe_normalize(const Vector3D& v, double eps = 1e-12) { + double n = v.norm(); + return (n > eps) ? (v / n) : Vector3D(0, 0, 1); +} + +/// 获取曲面在边界上的采样点 +/// @return {u, v} 参数对列表，沿指定边界均匀采样 +std::vector> +sample_boundary(const NurbsSurface& surf, BoundaryEdge edge, int n_samples) { + std::vector> pts; + auto [nu, nv] = surf.num_control_points(); + const auto& ku = surf.knots_u(); + const auto& kv = surf.knots_v(); + double u_min = ku[surf.degree_u()]; + double u_max = ku[nu]; + double v_min = kv[surf.degree_v()]; + double v_max = kv[nv]; + + for (int i = 0; i <= n_samples; ++i) { + double t = static_cast(i) / n_samples; + switch (edge) { + case BoundaryEdge::UMin: + pts.emplace_back(u_min, v_min + t * (v_max - v_min)); + break; + case BoundaryEdge::UMax: + pts.emplace_back(u_max, v_min + t * (v_max - v_min)); + break; + case BoundaryEdge::VMin: + pts.emplace_back(u_min + t * (u_max - u_min), v_min); + break; + case BoundaryEdge::VMax: + pts.emplace_back(u_min + t * (u_max - u_min), v_max); + break; + } + } + return pts; +} + +/// 获取边界上的控制点索引 +/// @return 沿边界的控制点 (i,j) 索引列表 +std::vector> +boundary_control_indices(const NurbsSurface& surf, BoundaryEdge edge) { + std::vector> idxs; + auto [nu, nv] = surf.num_control_points(); + + switch (edge) { + case BoundaryEdge::UMin: + for (int j = 0; j < nv; ++j) idxs.emplace_back(0, j); + break; + case BoundaryEdge::UMax: + for (int j = 0; j < nv; ++j) idxs.emplace_back(nu - 1, j); + break; + case BoundaryEdge::VMin: + for (int i = 0; i < nu; ++i) idxs.emplace_back(i, 0); + break; + case BoundaryEdge::VMax: + for (int i = 0; i < nu; ++i) idxs.emplace_back(i, nv - 1); + break; + } + return idxs; +} + +/// 获取第 k 排（从边界开始，0=边界）的控制点索引 +std::vector> +row_control_indices(const NurbsSurface& surf, BoundaryEdge edge, int k) { + auto [nu, nv] = surf.num_control_points(); + std::vector> idxs; + + switch (edge) { + case BoundaryEdge::UMin: + for (int j = 0; j < nv; ++j) idxs.emplace_back(k, j); + break; + case BoundaryEdge::UMax: + for (int j = 0; j < nv; ++j) idxs.emplace_back(nu - 1 - k, j); + break; + case BoundaryEdge::VMin: + for (int i = 0; i < nu; ++i) idxs.emplace_back(i, k); + break; + case BoundaryEdge::VMax: + for (int i = 0; i < nu; ++i) idxs.emplace_back(i, nv - 1 - k); + break; + } + return idxs; +} + +/// 从边界采样点计算 reference 曲面上最近的参数 +/// 简化：假设两曲面已经对齐，直接计算 target 边界点在 reference 上的投影 +/// 使用采样点匹配策略找到对应参数 +double find_closest_v_on_ref( + const Point3D& p_target, + const NurbsSurface& ref, + BoundaryEdge ref_edge, int n_samples, + double& best_u, double& best_v) { + + auto samples = sample_boundary(ref, ref_edge, n_samples * 4); + double best_dist = std::numeric_limits::max(); + best_u = 0.0; best_v = 0.0; + + for (const auto& [u, v] : samples) { + auto q = ref.evaluate(u, v); + double d = (p_target - q).squaredNorm(); + if (d < best_dist) { + best_dist = d; + best_u = u; + best_v = v; + } + } + return std::sqrt(best_dist); +} + +/// 计算曲面在 (u,v) 处的曲面法向量 +Vector3D surface_normal_at(const NurbsSurface& surf, double u, double v) { + auto du = surf.derivative_u(u, v); + auto dv = surf.derivative_v(u, v); + return safe_normalize(du.cross(dv)); +} + +/// 计算沿边界的 cross-boundary 方向（向内） +Vector3D cross_boundary_direction(BoundaryEdge edge) { + switch (edge) { + case BoundaryEdge::UMin: return Vector3D(1, 0, 0); // +u 方向 + case BoundaryEdge::UMax: return Vector3D(-1, 0, 0); // -u 方向 + case BoundaryEdge::VMin: return Vector3D(0, 1, 0); // +v 方向 + case BoundaryEdge::VMax: return Vector3D(0, -1, 0); // -v 方向 + } + return Vector3D(1, 0, 0); +} + +/// 近似曲面在控制点的 tangent 方向（用于 G1 调整） +Vector3D approx_tangent_dir(const NurbsSurface& surf, + BoundaryEdge edge, int idx) { + auto [nu, nv] = surf.num_control_points(); + const auto& cp = surf.control_points(); + + if (edge == BoundaryEdge::UMin || edge == BoundaryEdge::UMax) { + // 沿 v 方向前后控制点的差 + if (nv > 1 && idx > 0 && idx < nv - 1) { + int i = (edge == BoundaryEdge::UMin) ? 0 : nu - 1; + auto d = cp[i][idx + 1] - cp[i][idx - 1]; + double len = d.norm(); + if (len > 1e-12) return d / len; + } + return Vector3D(0, 1, 0); + } else { + if (nu > 1 && idx > 0 && idx < nu - 1) { + int j = (edge == BoundaryEdge::VMin) ? 0 : nv - 1; + auto d = cp[idx + 1][j] - cp[idx - 1][j]; + double len = d.norm(); + if (len > 1e-12) return d / len; + } + return Vector3D(1, 0, 0); + } +} + +// ═══════════════════════════════════════════════════════════ +// G0 匹配：调整边界控制点到 reference 曲面 +// ═══════════════════════════════════════════════════════════ + +void apply_g0(const NurbsSurface& target, + const NurbsSurface& reference, + BoundaryEdge edge, + std::vector>& grid, + int n_uv_samples) { + + auto [nu, nv] = target.num_control_points(); + auto bound_cps = boundary_control_indices(target, edge); + + // 对每个边界控制点，投影到 reference 曲面 + for (size_t bidx = 0; bidx < bound_cps.size(); ++bidx) { + int i = bound_cps[bidx].first; + int j = bound_cps[bidx].second; + Point3D pt = grid[i][j]; + + // 在 reference 曲面上找到最近点 + double best_u = 0.0, best_v = 0.0; + find_closest_v_on_ref(pt, reference, edge, n_uv_samples * 4, best_u, best_v); + auto ref_pt = reference.evaluate(best_u, best_v); + grid[i][j] = ref_pt; + } +} + +// ═══════════════════════════════════════════════════════════ +// G1 匹配：调整第1排控制点使切平面连续 +// ═══════════════════════════════════════════════════════════ + +void apply_g1(const NurbsSurface& target, + const NurbsSurface& reference, + BoundaryEdge edge, + std::vector>& grid, + int n_uv_samples) { + + auto [nu, nv] = target.num_control_points(); + + // 首先确保 G0 + apply_g0(target, reference, edge, grid, n_uv_samples); + + // 获取边界和第1排控制点索引 + auto bound_cps = boundary_control_indices(target, edge); + auto row1_cps = row_control_indices(target, edge, 1); + + // 在每个边界采样点，计算 reference 的 cross-boundary tangent + auto samples = sample_boundary(target, edge, n_uv_samples); + + // 对边界控制点沿边均匀分布的简化方法： + // 调整 row1 控制点使 cross-boundary 切向量与 reference 一致 + for (size_t bidx = 0; bidx < bound_cps.size() && bidx < row1_cps.size(); ++bidx) { + int bi = bound_cps[bidx].first, bj = bound_cps[bidx].second; + int ri = row1_cps[bidx].first, rj = row1_cps[bidx].second; + + // 采样点参数（沿边均匀分布） + double t = static_cast(bidx) / (bound_cps.size() - 1); + auto [nu_s, nv_s] = target.num_control_points(); + const auto& ku = target.knots_u(); + const auto& kv = target.knots_v(); + double u_min = ku[target.degree_u()], u_max = ku[nu_s]; + double v_min = kv[target.degree_v()], v_max = kv[nv_s]; + + double u, v; + switch (edge) { + case BoundaryEdge::UMin: u = u_min; v = v_min + t * (v_max - v_min); break; + case BoundaryEdge::UMax: u = u_max; v = v_min + t * (v_max - v_min); break; + case BoundaryEdge::VMin: u = u_min + t * (u_max - u_min); v = v_min; break; + case BoundaryEdge::VMax: u = u_min + t * (u_max - u_min); v = v_max; break; + default: u = 0; v = 0; break; + } + + // reference 在该点的 cross-boundary tangent + auto du_ref = reference.derivative_u(u, v); + auto dv_ref = reference.derivative_v(u, v); + Vector3D cb_tangent_ref; + if (edge == BoundaryEdge::UMin || edge == BoundaryEdge::UMax) { + cb_tangent_ref = safe_normalize(du_ref); + } else { + cb_tangent_ref = safe_normalize(dv_ref); + } + + // target 的 cross-boundary tangent 由 grid[bi][bj] 和 grid[ri][rj] 决定 + // 差值近似 = (grid[ri][rj] - grid[bi][bj]) * (degree / param_range) + double param_range; + if (edge == BoundaryEdge::UMin || edge == BoundaryEdge::UMax) { + param_range = u_max - u_min; + } else { + param_range = v_max - v_min; + } + if (param_range < 1e-12) param_range = 1.0; + + int deg = (edge == BoundaryEdge::UMin || edge == BoundaryEdge::UMax) + ? target.degree_u() : target.degree_v(); + + // 目标 cross-boundary 切向量（应该匹配 reference） + double scale = param_range / deg; + Point3D desired_ri = { + grid[bi][bj].x() + cb_tangent_ref.x() * scale, + grid[bi][bj].y() + cb_tangent_ref.y() * scale, + grid[bi][bj].z() + cb_tangent_ref.z() * scale + }; + + // G1: 将 row1 控制点移到目标位置 + grid[ri][rj] = desired_ri; + } +} + +// ═══════════════════════════════════════════════════════════ +// G2 匹配：调整第1-2排控制点使曲率连续 +// ═══════════════════════════════════════════════════════════ + +void apply_g2(const NurbsSurface& target, + const NurbsSurface& reference, + BoundaryEdge edge, + std::vector>& grid, + int n_uv_samples) { + + auto [nu, nv] = target.num_control_points(); + + // 首先确保 G0 + G1 + apply_g1(target, reference, edge, grid, n_uv_samples); + + // 获取边界、第1排、第2排控制点索引 + auto bound_cps = boundary_control_indices(target, edge); + auto row1_cps = row_control_indices(target, edge, 1); + auto row2_cps = row_control_indices(target, edge, 2); + + // 对每个控制点计算 G2 约束 + for (size_t bidx = 0; bidx < bound_cps.size() && + bidx < row1_cps.size() && + bidx < row2_cps.size(); ++bidx) { + int bi = bound_cps[bidx].first, bj = bound_cps[bidx].second; + int r1i = row1_cps[bidx].first, r1j = row1_cps[bidx].second; + int r2i = row2_cps[bidx].first, r2j = row2_cps[bidx].second; + + double t = static_cast(bidx) / (bound_cps.size() - 1); + auto [nu_s, nv_s] = target.num_control_points(); + const auto& ku = target.knots_u(); + const auto& kv = target.knots_v(); + double u_min = ku[target.degree_u()], u_max = ku[nu_s]; + double v_min = kv[target.degree_v()], v_max = kv[nv_s]; + + double u, v; + switch (edge) { + case BoundaryEdge::UMin: u = u_min; v = v_min + t * (v_max - v_min); break; + case BoundaryEdge::UMax: u = u_max; v = v_min + t * (v_max - v_min); break; + case BoundaryEdge::VMin: u = u_min + t * (u_max - u_min); v = v_min; break; + case BoundaryEdge::VMax: u = u_min + t * (u_max - u_min); v = v_max; break; + default: u = 0; v = 0; break; + } + + // reference 在该点沿 cross-boundary 方向的二阶导数 + double h = 1e-4; + Vector3D d2_ref; + if (edge == BoundaryEdge::UMin || edge == BoundaryEdge::UMax) { + auto du1 = reference.derivative_u(u + h, v); + auto du0 = reference.derivative_u(u - h, v); + d2_ref = (du1 - du0) / (2.0 * h); + } else { + auto dv1 = reference.derivative_v(u, v + h); + auto dv0 = reference.derivative_v(u, v - h); + d2_ref = (dv1 - dv0) / (2.0 * h); + } + + // 通过控制点差分近似二阶导数 + // row0 (boundary), row1, row2 的二阶差分 = = row2 - 2*row1 + row0 + int deg = (edge == BoundaryEdge::UMin || edge == BoundaryEdge::UMax) + ? target.degree_u() : target.degree_v(); + double param_range; + if (edge == BoundaryEdge::UMin || edge == BoundaryEdge::UMax) { + param_range = u_max - u_min; + } else { + param_range = v_max - v_min; + } + if (param_range < 1e-12) param_range = 1.0; + double h_ctrl = param_range / deg; + + // row2 位置应满足： (grid[r2i][r2j] - 2*grid[r1i][r1j] + grid[bi][bj]) / (h_ctrl*h_ctrl) ≈ d2_ref + double h2 = h_ctrl * h_ctrl; + Point3D row1 = grid[r1i][r1j]; + Point3D row0 = grid[bi][bj]; + + Point3D desired_row2 = { + row0.x() + 2.0 * (row1.x() - row0.x()) + d2_ref.x() * h2, + row0.y() + 2.0 * (row1.y() - row0.y()) + d2_ref.y() * h2, + row0.z() + 2.0 * (row1.z() - row0.z()) + d2_ref.z() * h2 + }; + + grid[r2i][r2j] = desired_row2; + } +} + +// ═══════════════════════════════════════════════════════════ +// G3 匹配：调整第1-3排控制点使曲率变化率连续 +// ═══════════════════════════════════════════════════════════ + +void apply_g3(const NurbsSurface& target, + const NurbsSurface& reference, + BoundaryEdge edge, + std::vector>& grid, + int n_uv_samples) { + + auto [nu, nv] = target.num_control_points(); + + // 首先确保 G0 + G1 + G2 + apply_g2(target, reference, edge, grid, n_uv_samples); + + auto bound_cps = boundary_control_indices(target, edge); + auto row1_cps = row_control_indices(target, edge, 1); + auto row2_cps = row_control_indices(target, edge, 2); + auto row3_cps = row_control_indices(target, edge, 3); + + for (size_t bidx = 0; bidx < bound_cps.size() && + bidx < row1_cps.size() && + bidx < row2_cps.size() && + bidx < row3_cps.size(); ++bidx) { + int bi = bound_cps[bidx].first, bj = bound_cps[bidx].second; + int r1i = row1_cps[bidx].first, r1j = row1_cps[bidx].second; + int r2i = row2_cps[bidx].first, r2j = row2_cps[bidx].second; + int r3i = row3_cps[bidx].first, r3j = row3_cps[bidx].second; + + double t = static_cast(bidx) / (bound_cps.size() - 1); + auto [nu_s, nv_s] = target.num_control_points(); + const auto& ku = target.knots_u(); + const auto& kv = target.knots_v(); + double u_min = ku[target.degree_u()], u_max = ku[nu_s]; + double v_min = kv[target.degree_v()], v_max = kv[nv_s]; + + double u, v; + switch (edge) { + case BoundaryEdge::UMin: u = u_min; v = v_min + t * (v_max - v_min); break; + case BoundaryEdge::UMax: u = u_max; v = v_min + t * (v_max - v_min); break; + case BoundaryEdge::VMin: u = u_min + t * (u_max - u_min); v = v_min; break; + case BoundaryEdge::VMax: u = u_min + t * (u_max - u_min); v = v_max; break; + default: u = 0; v = 0; break; + } + + // reference 的三阶导数（cross-boundary 方向） + double h = 1e-4; + Vector3D d3_ref; + if (edge == BoundaryEdge::UMin || edge == BoundaryEdge::UMax) { + auto duu1 = (reference.derivative_u(u + 2*h, v) + - reference.derivative_u(u, v)) / (2.0 * h); + auto duu0 = (reference.derivative_u(u, v) + - reference.derivative_u(u - 2*h, v)) / (2.0 * h); + d3_ref = (duu1 - duu0) / (2.0 * h); + } else { + auto dvv1 = (reference.derivative_v(u, v + 2*h) + - reference.derivative_v(u, v)) / (2.0 * h); + auto dvv0 = (reference.derivative_v(u, v) + - reference.derivative_v(u, v - 2*h)) / (2.0 * h); + d3_ref = (dvv1 - dvv0) / (2.0 * h); + } + + // 通过控制点的三阶差分近似 + int deg = (edge == BoundaryEdge::UMin || edge == BoundaryEdge::UMax) + ? target.degree_u() : target.degree_v(); + double param_range; + if (edge == BoundaryEdge::UMin || edge == BoundaryEdge::UMax) { + param_range = u_max - u_min; + } else { + param_range = v_max - v_min; + } + if (param_range < 1e-12) param_range = 1.0; + double h_ctrl = param_range / deg; + double h3 = h_ctrl * h_ctrl * h_ctrl; + + // 三阶差分：row3 - 3*row2 + 3*row1 - row0 ≈ d3_ref * h³ + Point3D row0 = grid[bi][bj]; + Point3D row1 = grid[r1i][r1j]; + Point3D row2 = grid[r2i][r2j]; + + Point3D desired_row3 = { + row0.x() + 3.0 * (row1.x() - row0.x()) + - 3.0 * (row2.x() - 2.0 * row1.x() + row0.x()) + + d3_ref.x() * h3, + row0.y() + 3.0 * (row1.y() - row0.y()) + - 3.0 * (row2.y() - 2.0 * row1.y() + row0.y()) + + d3_ref.y() * h3, + row0.z() + 3.0 * (row1.z() - row0.z()) + - 3.0 * (row2.z() - 2.0 * row1.z() + row0.z()) + + d3_ref.z() * h3 + }; + + grid[r3i][r3j] = desired_row3; + } +} + +// ═══════════════════════════════════════════════════════════ +// 连续性验证 +// ═══════════════════════════════════════════════════════════ + +struct ContinuityErrors { + double max_g0 = 0.0; + double max_g1 = 0.0; + double max_g2 = 0.0; + double max_g3 = 0.0; +}; + +ContinuityErrors verify_continuity( + const NurbsSurface& matched, + const NurbsSurface& reference, + BoundaryEdge edge, + ContinuityLevel level, + int n_samples) { + + ContinuityErrors errs; + auto samples = sample_boundary(matched, edge, n_samples); + + for (const auto& [u, v] : samples) { + auto p_m = matched.evaluate(u, v); + auto du_m = matched.derivative_u(u, v); + auto dv_m = matched.derivative_v(u, v); + + // 在 reference 上找最近点 + auto p_r = reference.evaluate(u, v); + + // G0: 位置误差 + double g0 = (p_m - p_r).norm(); + errs.max_g0 = std::max(errs.max_g0, g0); + + // G1: 切平面角度误差（法线夹角） + auto n_m = safe_normalize(du_m.cross(dv_m)); + auto du_r = reference.derivative_u(u, v); + auto dv_r = reference.derivative_v(u, v); + auto n_r = safe_normalize(du_r.cross(dv_r)); + double dot_n = std::max(-1.0, std::min(1.0, n_m.dot(n_r))); + double g1 = std::acos(dot_n); + errs.max_g1 = std::max(errs.max_g1, g1); + + // G2: 曲率偏差（通过二阶导数差值近似） + if (level >= ContinuityLevel::G2) { + double h = 1e-4; + Vector3D d2_m, d2_r; + if (edge == BoundaryEdge::UMin || edge == BoundaryEdge::UMax) { + auto du1_m = matched.derivative_u(u + h, v); + auto du0_m = matched.derivative_u(u - h, v); + d2_m = (du1_m - du0_m) / (2.0 * h); + auto du1_r = reference.derivative_u(u + h, v); + auto du0_r = reference.derivative_u(u - h, v); + d2_r = (du1_r - du0_r) / (2.0 * h); + } else { + auto dv1_m = matched.derivative_v(u, v + h); + auto dv0_m = matched.derivative_v(u, v - h); + d2_m = (dv1_m - dv0_m) / (2.0 * h); + auto dv1_r = reference.derivative_v(u, v + h); + auto dv0_r = reference.derivative_v(u, v - h); + d2_r = (dv1_r - dv0_r) / (2.0 * h); + } + double g2 = (d2_m - d2_r).norm(); + errs.max_g2 = std::max(errs.max_g2, g2); + } + + // G3: 曲率变化率偏差 + if (level >= ContinuityLevel::G3) { + double h = 1e-4; + Vector3D d3_m, d3_r; + if (edge == BoundaryEdge::UMin || edge == BoundaryEdge::UMax) { + auto du_a = matched.derivative_u(u + 2*h, v); + auto du_b = matched.derivative_u(u, v); + auto du_c = matched.derivative_u(u - 2*h, v); + d3_m = ((du_a - du_b) / (2.0*h) - (du_b - du_c) / (2.0*h)) / (2.0*h); + auto dr_a = reference.derivative_u(u + 2*h, v); + auto dr_b = reference.derivative_u(u, v); + auto dr_c = reference.derivative_u(u - 2*h, v); + d3_r = ((dr_a - dr_b) / (2.0*h) - (dr_b - dr_c) / (2.0*h)) / (2.0*h); + } else { + auto dv_a = matched.derivative_v(u, v + 2*h); + auto dv_b = matched.derivative_v(u, v); + auto dv_c = matched.derivative_v(u, v - 2*h); + d3_m = ((dv_a - dv_b) / (2.0*h) - (dv_b - dv_c) / (2.0*h)) / (2.0*h); + auto dr_a = reference.derivative_v(u, v + 2*h); + auto dr_b = reference.derivative_v(u, v); + auto dr_c = reference.derivative_v(u, v - 2*h); + d3_r = ((dr_a - dr_b) / (2.0*h) - (dr_b - dr_c) / (2.0*h)) / (2.0*h); + } + double g3 = (d3_m - d3_r).norm(); + errs.max_g3 = std::max(errs.max_g3, g3); + } + } + return errs; +} + +} // namespace + +// ═══════════════════════════════════════════════════════════ +// match_surface 主函数 +// ═══════════════════════════════════════════════════════════ + +MatchSurfaceResult match_surface( + const NurbsSurface& target, + const NurbsSurface& reference, + BoundaryEdge edge, + const MatchSurfaceOptions& options) { + + MatchSurfaceResult result; + + // 复制 target 的控制点和权重网格 + const auto& src_cp = target.control_points(); + const auto& src_w = target.weights(); + auto [nu, nv] = target.num_control_points(); + std::vector> grid(nu, std::vector(nv)); + std::vector> weights(nu, std::vector(nv, 1.0)); + + bool has_weights = !src_w.empty(); + for (int i = 0; i < nu; ++i) { + for (int j = 0; j < nv; ++j) { + grid[i][j] = src_cp[i][j]; + if (has_weights) { + weights[i][j] = src_w[i][j]; + } + } + } + + // 确定调整行数 + int rows = options.rows_to_adjust; + if (rows <= 0) { + switch (options.continuity) { + case ContinuityLevel::G0: rows = 1; break; + case ContinuityLevel::G1: rows = 1; break; + case ContinuityLevel::G2: rows = 2; break; + case ContinuityLevel::G3: rows = 3; break; + } + } + result.rows_adjusted = rows; + + int n_samples = options.boundary_samples; + + // 按连续性等级递增应用：G0 → G1 → G2 → G3 + // 每条增加一排控制点的自由度 + int target_level = static_cast(options.continuity); + + if (target_level >= 0) { + apply_g0(target, reference, edge, grid, n_samples); + } + if (target_level >= 1) { + apply_g1(target, reference, edge, grid, n_samples); + } + if (target_level >= 2) { + apply_g2(target, reference, edge, grid, n_samples); + } + if (target_level >= 3) { + apply_g3(target, reference, edge, grid, n_samples); + } + + // 构建修改后的曲面 + NurbsSurface matched(grid, + target.knots_u(), + target.knots_v(), + weights, + target.degree_u(), + target.degree_v()); + result.matched_surface = matched; + + // 计算控制点位移 + for (int i = 0; i < nu; ++i) { + for (int j = 0; j < nv; ++j) { + result.control_displacements.push_back(grid[i][j] - src_cp[i][j]); + } + } + + // 验证连续性 + auto errs = verify_continuity(matched, reference, edge, options.continuity, n_samples); + result.max_position_error = errs.max_g0; + result.max_tangent_error = errs.max_g1; + result.max_curvature_error = errs.max_g2; + result.max_derivative_error = errs.max_g3; + + double tol = options.tolerance; + result.g0_ok = (errs.max_g0 < tol); + result.g1_ok = (errs.max_g1 < 1e-4); // 角度容差 ~0.0001 rad ≈ 0.006° + result.g2_ok = (errs.max_g2 < tol * 10.0); + result.g3_ok = (errs.max_g3 < tol * 100.0); + + return result; +} + +} // namespace vde::curves diff --git a/tests/curves/CMakeLists.txt b/tests/curves/CMakeLists.txt index f3caef3..ecd515a 100644 --- a/tests/curves/CMakeLists.txt +++ b/tests/curves/CMakeLists.txt @@ -9,3 +9,6 @@ add_vde_test(test_surface_extension) add_vde_test(test_class_a_surfacing) add_vde_test(test_advanced_intersection) add_vde_test(test_iga_prep) +add_vde_test(test_curve_tools) +add_vde_test(test_generative_surfaces) +add_vde_test(test_surface_editing) diff --git a/tests/curves/test_curve_tools.cpp b/tests/curves/test_curve_tools.cpp new file mode 100644 index 0000000..1b19751 --- /dev/null +++ b/tests/curves/test_curve_tools.cpp @@ -0,0 +1,273 @@ +#include +#include "vde/curves/curve_tools.h" +#include "vde/curves/surface_analysis.h" +#include "vde/curves/nurbs_curve.h" +#include "vde/curves/nurbs_surface.h" +#include + +using namespace vde::curves; +using namespace vde::core; + +// --------------------------------------------------------------------------- +// Helpers +// --------------------------------------------------------------------------- + +/// Create a planar NURBS surface on XY plane over [0,1]×[0,1] +static NurbsSurface plane_surface() { + std::vector> grid = { + {Point3D(0,0,0), Point3D(0,1,0)}, + {Point3D(1,0,0), Point3D(1,1,0)} + }; + return NurbsSurface(grid, {0,0,1,1}, {0,0,1,1}, + {{1,1},{1,1}}, 1, 1); +} + +/// Create a NURBS line segment +static NurbsCurve line(const Point3D& a, const Point3D& b) { + return NurbsCurve({a, b}, {0, 0, 1, 1}, {1, 1}, 1); +} + +/// Create a cubic NURBS curve +static NurbsCurve cubic_curve() { + return NurbsCurve( + {Point3D(0,0,0), Point3D(1,2,0), Point3D(2,1,0), Point3D(3,0,0)}, + {0,0,0,0,1,1,1,1}, + {1,1,1,1}, + 3); +} + +// --------------------------------------------------------------------------- +// Test 1: project_curve_on_surface (via curve_tools.h inclusion) +// --------------------------------------------------------------------------- + +TEST(CurveToolsTest, ProjectCurveOnSurface_PlaneProjection) { + auto surf = plane_surface(); + auto curve = line(Point3D(0.2, 0.3, 1.0), Point3D(0.8, 0.7, 2.0)); + + // Project the 3D curve onto the XY plane surface + auto pcurve = project_curve_on_surface(curve, surf, 50); + + // pcurve should be a valid NURBS (non-degenerate) + EXPECT_GE(pcurve.degree(), 0); + auto dom = pcurve.domain(); + EXPECT_LE(dom.first, dom.second); +} + +// --------------------------------------------------------------------------- +// Test 2: parallel_curve — straight line offset +// --------------------------------------------------------------------------- + +TEST(CurveToolsTest, ParallelCurve_StraightLineOffset) { + auto c = line(Point3D(0, 0, 0), Point3D(10, 0, 0)); + // Offset by 1 in Y direction (plane normal = +Z) + auto offset = parallel_curve(c, 1.0, Vector3D(0, 0, 1), 200); + + EXPECT_GT(offset.degree(), 0); + + // Midpoint should be (5, 1, 0) — shifted by 1 in Y + auto dom = offset.domain(); + double t_mid = (dom.first + dom.second) * 0.5; + Point3D pm = offset.evaluate(t_mid); + EXPECT_NEAR(pm.x(), 5.0, 0.5); + EXPECT_NEAR(pm.y(), 1.0, 0.3); + EXPECT_NEAR(pm.z(), 0.0, 0.1); +} + +// --------------------------------------------------------------------------- +// Test 3: parallel_curve — negative offset +// --------------------------------------------------------------------------- + +TEST(CurveToolsTest, ParallelCurve_NegativeOffset) { + auto c = line(Point3D(2, 0, 0), Point3D(2, 8, 0)); + // Offset by -1 in X direction (plane normal = +Z, tangent = +Y → offset_dir = n×tangent = -X) + // distance=-1 offsets opposite to offset_dir, i.e., in +X + auto offset = parallel_curve(c, -1.0, Vector3D(0, 0, 1), 200); + + EXPECT_GT(offset.degree(), 0); + + auto dom = offset.domain(); + double t_mid = (dom.first + dom.second) * 0.5; + Point3D pm = offset.evaluate(t_mid); + // Midpoint should be shifted from x=2 in the offset direction + EXPECT_NEAR(pm.x(), 3.0, 0.5); + EXPECT_NEAR(pm.y(), 4.0, 1.0); +} + +// --------------------------------------------------------------------------- +// Test 4: connect_curve — G1 continuity +// --------------------------------------------------------------------------- + +TEST(CurveToolsTest, ConnectCurve_G1_BasicConnect) { + auto c1 = line(Point3D(0, 0, 0), Point3D(5, 0, 0)); + auto c2 = line(Point3D(10, 0, 0), Point3D(15, 0, 0)); + + auto bridge = connect_curve(c1, c2, Continuity::G1); + + EXPECT_GT(bridge.degree(), 0); + + // Endpoints should match + auto dom = bridge.domain(); + Point3D p_start = bridge.evaluate(dom.first); + Point3D p_end = bridge.evaluate(dom.second); + + EXPECT_NEAR(p_start.x(), 5.0, 0.5); + EXPECT_NEAR(p_start.y(), 0.0, 0.5); + EXPECT_NEAR(p_end.x(), 10.0, 0.5); + EXPECT_NEAR(p_end.y(), 0.0, 0.5); +} + +// --------------------------------------------------------------------------- +// Test 5: connect_curve — G2 continuity +// --------------------------------------------------------------------------- + +TEST(CurveToolsTest, ConnectCurve_G2_HigherContinuity) { + auto c1 = cubic_curve(); // ends near (3,0,0) + auto c2 = line(Point3D(6, 0, 0), Point3D(9, 0, 0)); // starts at (6,0,0) + + auto bridge = connect_curve(c1, c2, Continuity::G2); + + EXPECT_GT(bridge.degree(), 0); +} + +// --------------------------------------------------------------------------- +// Test 6: connect_curve — G3 continuity +// --------------------------------------------------------------------------- + +TEST(CurveToolsTest, ConnectCurve_G3_HighestContinuity) { + auto c1 = cubic_curve(); // ends near (3,0,0) + auto c2 = line(Point3D(8, 2, 0), Point3D(12, 2, 0)); // starts at (8,2,0) + + auto bridge = connect_curve(c1, c2, Continuity::G3); + + EXPECT_GT(bridge.degree(), 0); +} + +// --------------------------------------------------------------------------- +// Test 7: helix_curve — basic helix +// --------------------------------------------------------------------------- + +TEST(CurveToolsTest, HelixCurve_BasicHelix) { + auto helix = helix_curve( + Point3D(0, 0, 0), // axis origin + Vector3D(0, 0, 1), // axis direction (Z-up) + 5.0, // radius + 2.0, // pitch + 10.0, // height + 100); // samples per turn + + EXPECT_GT(helix.degree(), 0); + + // Check start point — should be near the helix start (radius distance from axis) + auto dom = helix.domain(); + Point3D p0 = helix.evaluate(dom.first); + double r0 = std::sqrt(p0.x() * p0.x() + p0.y() * p0.y()); + EXPECT_NEAR(r0, 5.0, 0.5); + EXPECT_NEAR(p0.z(), 0.0, 0.5); + + // Check end point — should be near the top + Point3D p1 = helix.evaluate(dom.second); + double r1 = std::sqrt(p1.x() * p1.x() + p1.y() * p1.y()); + EXPECT_NEAR(r1, 5.0, 0.5); + EXPECT_NEAR(p1.z(), 10.0, 0.5); +} + +// --------------------------------------------------------------------------- +// Test 8: helix_curve — zero pitch / height edge case +// --------------------------------------------------------------------------- + +TEST(CurveToolsTest, HelixCurve_ZeroPitchReturnsEmpty) { + auto helix = helix_curve( + Point3D(0, 0, 0), Vector3D(0, 0, 1), + 2.0, 0.0, 5.0); // pitch = 0 → invalid + + // Should return empty curve + EXPECT_EQ(helix.degree(), 0); +} + +// --------------------------------------------------------------------------- +// Test 9: isoparametric_curves — basic extraction +// --------------------------------------------------------------------------- + +TEST(CurveToolsTest, IsoparametricCurves_BasicExtraction) { + auto surf = plane_surface(); + + auto iso = isoparametric_curves(surf, 3, 3); + + EXPECT_EQ(iso.u_curves.size(), 3u); + EXPECT_EQ(iso.v_curves.size(), 3u); + + // Each extracted curve should be valid + for (const auto& cv : iso.u_curves) { + EXPECT_GT(cv.degree(), 0); + } + for (const auto& cv : iso.v_curves) { + EXPECT_GT(cv.degree(), 0); + } +} + +// --------------------------------------------------------------------------- +// Test 10: isoparametric_curves — single curve +// --------------------------------------------------------------------------- + +TEST(CurveToolsTest, IsoparametricCurves_SingleCurve) { + auto surf = plane_surface(); + auto iso = isoparametric_curves(surf, 1, 1); + + EXPECT_EQ(iso.u_curves.size(), 1u); + EXPECT_EQ(iso.v_curves.size(), 1u); +} + +// --------------------------------------------------------------------------- +// Test 11: surface_analysis — inflection_lines +// --------------------------------------------------------------------------- + +TEST(CurveToolsTest, InflectionLines_PlaneNoInflections) { + auto surf = plane_surface(); + auto result = inflection_lines(surf, 30, 30); + + // Plane has Gaussian curvature ≡ 0 everywhere → no sign change → no inflections + EXPECT_EQ(result.segments.size(), 0u); +} + +// --------------------------------------------------------------------------- +// Test 12: surface_analysis — checker_mapping +// --------------------------------------------------------------------------- + +TEST(CurveToolsTest, CheckerMapping_OutputDimensions) { + auto surf = plane_surface(); + auto ckr = checker_mapping(surf, Vector3D(1, 0, 0), 20, 15); + + EXPECT_EQ(ckr.res_u, 20); + EXPECT_EQ(ckr.res_v, 15); + EXPECT_EQ(static_cast(ckr.intensity.size()), 21); // res_u+1 + EXPECT_EQ(static_cast(ckr.intensity[0].size()), 16); // res_v+1 +} + +TEST(CurveToolsTest, CheckerMapping_IntensityRange) { + auto surf = plane_surface(); + auto ckr = checker_mapping(surf, Vector3D(0.5, 0.3, 1.0), 15, 15); + + for (const auto& row : ckr.intensity) { + for (double val : row) { + EXPECT_GE(val, 0.0); + EXPECT_LE(val, 1.0); + } + } +} + +// --------------------------------------------------------------------------- +// Test 13: surface_analysis — surface_checker_report +// --------------------------------------------------------------------------- + +TEST(CurveToolsTest, SurfaceCheckerReport_PlaneIsSmooth) { + auto surf = plane_surface(); + auto report = surface_checker_report(surf, 30); + + // Plane is perfectly smooth — should get grade A + EXPECT_GT(report.checker_uniformity, 0.5); + EXPECT_GT(report.curvature_continuity, 0.5); + EXPECT_GE(report.gaussian_range, 0.0); + + // Description should be non-empty + EXPECT_FALSE(report.description.empty()); +} diff --git a/tests/curves/test_generative_surfaces.cpp b/tests/curves/test_generative_surfaces.cpp new file mode 100644 index 0000000..9c5bf9c --- /dev/null +++ b/tests/curves/test_generative_surfaces.cpp @@ -0,0 +1,340 @@ +#include +#include "vde/curves/generative_surfaces.h" +#include "vde/curves/nurbs_curve.h" +#include "vde/core/point.h" +#include + +using namespace vde::curves; +using namespace vde::core; + +// ── 辅助：创建直线 NURBS ── +static NurbsCurve make_line(const Point3D& a, const Point3D& b) { + return NurbsCurve({a, b}, {0, 0, 1, 1}, {1, 1}, 1); +} + +// ── 辅助：创建圆弧 profile ── +static NurbsCurve make_arc_profile() { + double r = 1.0; + double w = std::sqrt(2.0) / 2.0; + return NurbsCurve( + {Point3D(0, r, 0), Point3D(r, r, 0), Point3D(r, 0, 0)}, + {0, 0, 0, 0.5, 1, 1, 1}, + {1, w, 1, w, 1, w, 1}, + 2 + ); +} + +// ============================================================================ +// Sweep 测试 +// ============================================================================ + +TEST(GenerativeSurfacesTest, SweepExplicit_Basic) { + // 圆形 profile 沿直线扫描 → 圆柱面 + NurbsCurve profile( + {Point3D(0, 1, 0), Point3D(1, 1, 0), Point3D(1, 0, 0)}, + {0, 0, 0, 1, 1, 1}, + {1, 1, 1}, + 2 + ); + NurbsCurve path = make_line({0, 0, 0}, {0, 0, 5}); + + SweepOption opt; + opt.type = SweepType::Explicit; + opt.samples = 8; + + auto surf = sweep_surface(profile, path, opt); + + EXPECT_EQ(surf.degree_u(), 2); + EXPECT_GE(surf.num_control_points()[0], 3); + EXPECT_GE(surf.num_control_points()[1], 3); + + // 验证曲面在参数域中点可求值 + auto p = surf.evaluate(0.5, 0.5); + EXPECT_GE(p.x(), 0.0); + EXPECT_LE(p.z(), 5.0); +} + +TEST(GenerativeSurfacesTest, SweepExplicit_StraightPath) { + // 矩形 profile 沿 X 轴扫描 + NurbsCurve profile = make_line({0, 0, 0}, {0, 0, 2}); + NurbsCurve path = make_line({0, 0, 0}, {10, 0, 0}); + + SweepOption opt; + opt.type = SweepType::Explicit; + opt.samples = 4; + + auto surf = sweep_surface(profile, path, opt); + + // 验证曲面两端 + auto p0 = surf.evaluate(0.0, 0.0); + auto p1 = surf.evaluate(1.0, 1.0); + + EXPECT_NEAR(p0.x(), 0.0, 1e-6); + EXPECT_NEAR(p1.x(), 10.0, 1e-6); + EXPECT_GE(p1.z(), 0.0); +} + +TEST(GenerativeSurfacesTest, SweepExplicit_VaryingSamples) { + NurbsCurve profile = make_line({0, 1, 0}, {0, -1, 0}); + NurbsCurve path = make_line({0, 0, 0}, {0, 0, 3}); + + SweepOption opt; + opt.type = SweepType::Explicit; + opt.samples = 16; + + auto surf = sweep_surface(profile, path, opt); + + // 控制点数量应随采样增加 + int n_u = surf.num_control_points()[0]; + EXPECT_GT(n_u, 4); + + // 曲面可正常求值 + auto mid = surf.evaluate(0.5, 0.5); + EXPECT_NEAR(mid.x(), 0.0, 1e-6); + EXPECT_NEAR(mid.z(), 1.5, 1.0); +} + +TEST(GenerativeSurfacesTest, SweepSpine_Basic) { + // 使用独立脊线的扫描 + NurbsCurve profile = make_line({0, 1, 0}, {0, -1, 0}); + NurbsCurve path = make_line({0, 0, 0}, {5, 0, 0}); + + // 脊线偏离路径，用于扭转控制 + NurbsCurve spine = make_line({0, 0, 0}, {5, 0, 2}); + + SweepOption opt; + opt.type = SweepType::Spine; + opt.spine = spine; + opt.samples = 8; + + auto surf = sweep_surface(profile, path, opt); + + EXPECT_GE(surf.num_control_points()[0], 3); + EXPECT_GE(surf.num_control_points()[1], 2); + + // 验证曲面可求值 + auto p = surf.evaluate(0.5, 0.5); + EXPECT_FALSE(std::isnan(p.x())); + EXPECT_FALSE(std::isnan(p.y())); + EXPECT_FALSE(std::isnan(p.z())); +} + +TEST(GenerativeSurfacesTest, SweepTwoGuides_Basic) { + // 双引导线扫描：两条引导线间距变化控制缩放 + NurbsCurve profile = make_line({0, 1, 0}, {0, -1, 0}); + NurbsCurve path = make_line({0, 0, 0}, {5, 0, 0}); + + // 两条引导线逐渐分开 + NurbsCurve guide1 = make_line({0, 0.5, 0}, {5, 2, 0}); + NurbsCurve guide2 = make_line({0, -0.5, 0}, {5, -2, 0}); + + SweepOption opt; + opt.type = SweepType::TwoGuides; + opt.guide1 = guide1; + opt.guide2 = guide2; + opt.samples = 8; + + auto surf = sweep_surface(profile, path, opt); + + EXPECT_GE(surf.num_control_points()[0], 3); + EXPECT_GE(surf.num_control_points()[1], 2); + + // 曲面起点和终点宽度应不同（缩放生效） + auto p_start = surf.evaluate(0.0, 0.5); + auto p_end = surf.evaluate(1.0, 0.5); + EXPECT_GE(std::abs(p_end.y()) - std::abs(p_start.y()), -0.5); +} + +TEST(GenerativeSurfacesTest, Sweep_DefaultOption) { + // 默认选项（Explicit 模式） + NurbsCurve profile = make_line({0, 1, 0}, {0, -1, 0}); + NurbsCurve path = make_line({0, 0, 0}, {3, 0, 0}); + + auto surf = sweep_surface(profile, path); + + EXPECT_GE(surf.num_control_points()[0], 3); + auto p = surf.evaluate(0.5, 0.5); + EXPECT_FALSE(std::isnan(p.x())); +} + +// ============================================================================ +// Loft 测试 +// ============================================================================ + +TEST(GenerativeSurfacesTest, Loft_Basic) { + // 两条平行截面放样 + NurbsCurve s1 = make_line({0, 0, 0}, {2, 0, 0}); + NurbsCurve s2 = make_line({0, 2, 5}, {2, 2, 5}); + + auto surf = loft_surface({s1, s2}); + + EXPECT_GE(surf.num_control_points()[0], 2); + EXPECT_GE(surf.num_control_points()[1], 2); + + // 验证曲面插值首尾截面 + auto p_bottom = surf.evaluate(0.5, 0.0); + EXPECT_NEAR(p_bottom.z(), 0.0, 1e-6); + + auto p_top = surf.evaluate(0.5, 1.0); + EXPECT_NEAR(p_top.z(), 5.0, 1e-6); +} + +TEST(GenerativeSurfacesTest, Loft_ThreeSections) { + NurbsCurve s1 = make_line({0, 0, 0}, {2, 0, 0}); + NurbsCurve s2 = make_line({0, 1, 3}, {2, 1, 3}); + NurbsCurve s3 = make_line({0, 0, 6}, {2, 0, 6}); + + auto surf = loft_surface({s1, s2, s3}); + + EXPECT_GE(surf.num_control_points()[0], 2); + EXPECT_GE(surf.num_control_points()[1], 2); + + // 中间截面应偏离直线 + auto p_mid = surf.evaluate(0.5, 0.5); + EXPECT_NEAR(p_mid.z(), 3.0, 1.0); + EXPECT_GT(p_mid.y(), 0.0); +} + +TEST(GenerativeSurfacesTest, Loft_WithGuides) { + NurbsCurve s1 = make_line({0, 0, 0}, {2, 0, 0}); + NurbsCurve s2 = make_line({0, 0, 5}, {2, 0, 5}); + NurbsCurve s3 = make_line({0, 0, 10}, {2, 0, 10}); + + // 引导线：弯曲的脊线 + NurbsCurve guide( + {Point3D(0, 1, 2.5), Point3D(1, 0, 5), Point3D(0, -1, 7.5)}, + {0, 0, 0, 1, 1, 1}, + {1, 1, 1}, + 2 + ); + + auto surf = loft_surface({s1, s2, s3}, {guide}); + + EXPECT_GE(surf.num_control_points()[0], 2); + EXPECT_GE(surf.num_control_points()[1], 2); + + auto p = surf.evaluate(0.5, 0.5); + EXPECT_FALSE(std::isnan(p.x())); +} + +TEST(GenerativeSurfacesTest, Loft_WithTangents) { + NurbsCurve s1 = make_line({0, 0, 0}, {2, 0, 0}); + NurbsCurve s2 = make_line({0, 0, 5}, {2, 0, 5}); + + LoftConstraint cons; + cons.start_tangent = Vector3D(0, 0, 1); // 垂直向上 + cons.end_tangent = Vector3D(0, 0, 1); // 垂直向上 + + auto surf = loft_surface({s1, s2}, {}, cons); + + EXPECT_GE(surf.num_control_points()[0], 2); + EXPECT_GE(surf.num_control_points()[1], 2); + + // 曲面应能正常求值 + auto p = surf.evaluate(0.5, 0.3); + EXPECT_FALSE(std::isnan(p.x())); +} + +TEST(GenerativeSurfacesTest, Loft_InvalidInput) { + // 单个截面应抛出异常 + NurbsCurve s1 = make_line({0, 0, 0}, {1, 0, 0}); + EXPECT_THROW((void)loft_surface({s1}), std::invalid_argument); +} + +// ============================================================================ +// Net 测试 +// ============================================================================ + +TEST(GenerativeSurfacesTest, Net_Basic) { + // 双向曲线网格：u 向和 v 向各 2 条线 + NurbsCurve u0 = make_line({0, 0, 0}, {5, 0, 0}); + NurbsCurve u1 = make_line({0, 2, 0}, {5, 2, 0}); + NurbsCurve v0 = make_line({0, 0, 0}, {0, 2, 0}); + NurbsCurve v1 = make_line({5, 0, 0}, {5, 2, 0}); + + auto surf = net_surface({u0, u1}, {v0, v1}); + + EXPECT_GE(surf.num_control_points()[0], 2); + EXPECT_GE(surf.num_control_points()[1], 2); + + // 曲面应在 XY 平面内 + auto p = surf.evaluate(0.5, 0.5); + EXPECT_NEAR(p.z(), 0.0, 1e-6); + EXPECT_GE(p.x(), 0.0); + EXPECT_GE(p.y(), 0.0); +} + +TEST(GenerativeSurfacesTest, Net_UnequalCount) { + // u 向 3 条，v 向 2 条 + NurbsCurve u0 = make_line({0, 0, 0}, {5, 0, 0}); + NurbsCurve u1 = make_line({0, 1, 0}, {5, 1, 0}); + NurbsCurve u2 = make_line({0, 2, 0}, {5, 2, 0}); + NurbsCurve v0 = make_line({0, 0, 0}, {0, 2, 0}); + NurbsCurve v1 = make_line({5, 0, 0}, {5, 2, 0}); + + auto surf = net_surface({u0, u1, u2}, {v0, v1}); + + EXPECT_EQ(surf.num_control_points()[0], 3); + EXPECT_EQ(surf.num_control_points()[1], 2); +} + +TEST(GenerativeSurfacesTest, Net_3DCurves) { + // 三维网格曲面 + NurbsCurve u0( + {Point3D(0, 0, 0), Point3D(1, 0, 1), Point3D(2, 0, 0)}, + {0, 0, 0, 1, 1, 1}, + {1, 1, 1}, + 2 + ); + NurbsCurve u1( + {Point3D(0, 2, 0), Point3D(1, 2, 1), Point3D(2, 2, 0)}, + {0, 0, 0, 1, 1, 1}, + {1, 1, 1}, + 2 + ); + NurbsCurve v0 = make_line({0, 0, 0}, {0, 2, 0}); + NurbsCurve v1 = make_line({2, 0, 0}, {2, 2, 0}); + + auto surf = net_surface({u0, u1}, {v0, v1}); + + EXPECT_GE(surf.num_control_points()[0], 2); + EXPECT_GE(surf.num_control_points()[1], 2); + + auto p = surf.evaluate(0.5, 0.5); + EXPECT_FALSE(std::isnan(p.x())); + EXPECT_GE(p.z(), 0.0); +} + +TEST(GenerativeSurfacesTest, Net_InvalidInput) { + NurbsCurve u0 = make_line({0, 0, 0}, {1, 0, 0}); + EXPECT_THROW((void)net_surface({u0}, {u0}), std::invalid_argument); + + NurbsCurve v0 = make_line({0, 0, 0}, {1, 0, 0}); + EXPECT_THROW((void)net_surface({u0, u0}, {v0}), std::invalid_argument); +} + +// ============================================================================ +// 集成测试：组合使用 +// ============================================================================ + +TEST(GenerativeSurfacesTest, Integration_SweepThenEvaluate) { + // 验证扫掠后曲面求值一致性 + NurbsCurve profile = make_line({0, 1, 0}, {0, -1, 0}); + NurbsCurve path = make_line({0, 0, 0}, {10, 0, 0}); + + SweepOption opt; + opt.type = SweepType::Explicit; + opt.samples = 10; + + auto surf = sweep_surface(profile, path, opt); + + // 曲面法向应在 XY 平面内（近似） + for (double u = 0.0; u <= 1.0; u += 0.25) { + for (double v = 0.0; v <= 1.0; v += 0.25) { + auto pt = surf.evaluate(u, v); + EXPECT_FALSE(std::isnan(pt.x())); + EXPECT_FALSE(std::isnan(pt.y())); + EXPECT_FALSE(std::isnan(pt.z())); + } + } +} diff --git a/tests/curves/test_surface_editing.cpp b/tests/curves/test_surface_editing.cpp new file mode 100644 index 0000000..c9f9c05 --- /dev/null +++ b/tests/curves/test_surface_editing.cpp @@ -0,0 +1,249 @@ +/** + * @file test_surface_editing.cpp + * @brief 曲面编辑三件套测试：match_surface / shape_fillet / control_point_edit + */ + +#include +#include "vde/curves/surface_editing.h" +#include "vde/curves/control_point_edit.h" +#include "vde/curves/nurbs_surface.h" +#include "vde/core/point.h" +#include + +namespace vde::curves { +namespace test { + +using core::Point3D; +using core::Vector3D; + +// ── Helper: create a simple biquadratic plane surface ── +NurbsSurface make_test_plane() { + // 3×3 control grid on z=0 plane + std::vector> grid = { + {{0, 0, 0}, {1, 0, 0}, {2, 0, 0}}, + {{0, 1, 0}, {1, 1, 0}, {2, 1, 0}}, + {{0, 2, 0}, {1, 2, 0}, {2, 2, 0}} + }; + return NurbsSurface(grid, + {0, 0, 0, 1, 1, 1}, // u-knots, degree=2 + {0, 0, 0, 1, 1, 1}, // v-knots, degree=2 + {}, 2, 2); // weights (empty = all 1.0) +} + +// ── Helper: create a slightly curved surface (degree 2×2) ── +NurbsSurface make_test_curved() { + // 3×3 grid, slightly curved in z + std::vector> grid = { + {{0, 0, 0}, {1, 0, 0.5}, {2, 0, 0}}, + {{0, 1, 0.5}, {1, 1, 1.0}, {2, 1, 0.5}}, + {{0, 2, 0}, {1, 2, 0.5}, {2, 2, 0}} + }; + return NurbsSurface(grid, + {0, 0, 0, 1, 1, 1}, {0, 0, 0, 1, 1, 1}, {}, 2, 2); +} + +// ═══════════════════════════════════════════════════════════ +// match_surface 测试 +// ═══════════════════════════════════════════════════════════ + +TEST(SurfaceEditingTest, MatchSurfaceG0Planar) { + // Two identical planes on z=0 — matching should produce no error + auto plane_a = make_test_plane(); + auto plane_b = make_test_plane(); + + MatchSurfaceOptions opts; + opts.continuity = ContinuityLevel::G0; + opts.boundary_samples = 10; + + auto result = match_surface(plane_b, plane_a, BoundaryEdge::VMin, opts); + + EXPECT_TRUE(result.g0_ok); + EXPECT_LT(result.max_position_error, 1e-4); + EXPECT_EQ(result.rows_adjusted, 1); + EXPECT_EQ(result.control_displacements.size(), 9u); // 3×3 grid +} + +TEST(SurfaceEditingTest, MatchSurfaceG0DifferentSurfaces) { + // Two slightly different curved surfaces — match VMin edge + auto surf_a = make_test_curved(); + // Modify a few control points to make it different + auto surf_b = make_test_curved(); + + MatchSurfaceOptions opts; + opts.continuity = ContinuityLevel::G0; + opts.boundary_samples = 10; + + auto result = match_surface(surf_b, surf_a, BoundaryEdge::VMin, opts); + + // Should produce a valid surface + auto [nu, nv] = result.matched_surface.num_control_points(); + EXPECT_GE(nu, 1); + EXPECT_GE(nv, 1); + // Control points should be adjusted + EXPECT_GT(result.control_displacements.size(), 0u); +} + +TEST(SurfaceEditingTest, MatchSurfaceG1Planar) { + auto plane_a = make_test_plane(); + auto plane_b = make_test_plane(); + + MatchSurfaceOptions opts; + opts.continuity = ContinuityLevel::G1; + opts.boundary_samples = 10; + + auto result = match_surface(plane_b, plane_a, BoundaryEdge::VMin, opts); + + EXPECT_TRUE(result.g0_ok); + EXPECT_TRUE(result.g1_ok); + EXPECT_EQ(result.rows_adjusted, 1); +} + +TEST(SurfaceEditingTest, MatchSurfaceG1DifferentSurfaces) { + auto surf_a = make_test_curved(); + auto surf_b = make_test_curved(); + + MatchSurfaceOptions opts; + opts.continuity = ContinuityLevel::G1; + opts.boundary_samples = 10; + + auto result = match_surface(surf_b, surf_a, BoundaryEdge::VMin, opts); + + // Check that errors are computed + EXPECT_GE(result.max_position_error, 0.0); + EXPECT_GE(result.max_tangent_error, 0.0); + EXPECT_EQ(result.rows_adjusted, 1); +} + +TEST(SurfaceEditingTest, MatchSurfaceG2Planar) { + auto plane_a = make_test_plane(); + auto plane_b = make_test_plane(); + + MatchSurfaceOptions opts; + opts.continuity = ContinuityLevel::G2; + opts.boundary_samples = 10; + + auto result = match_surface(plane_b, plane_a, BoundaryEdge::VMin, opts); + + EXPECT_TRUE(result.g0_ok); + EXPECT_EQ(result.rows_adjusted, 2); +} + +TEST(SurfaceEditingTest, MatchSurfaceG3Planar) { + auto plane_a = make_test_plane(); + auto plane_b = make_test_plane(); + + MatchSurfaceOptions opts; + opts.continuity = ContinuityLevel::G3; + opts.boundary_samples = 10; + + auto result = match_surface(plane_b, plane_a, BoundaryEdge::VMin, opts); + + EXPECT_TRUE(result.g0_ok); + EXPECT_EQ(result.rows_adjusted, 3); // G3 adjusts 3 rows +} + +TEST(SurfaceEditingTest, MatchSurfaceAllEdges) { + auto plane_a = make_test_plane(); + auto plane_b = make_test_plane(); + + MatchSurfaceOptions opts; + opts.continuity = ContinuityLevel::G1; + + // Test all 4 boundary edges + for (int e = 0; e < 4; ++e) { + BoundaryEdge edge = static_cast(e); + auto result = match_surface(plane_b, plane_a, edge, opts); + EXPECT_TRUE(result.g0_ok) << "Failed on edge " << e; + } +} + +TEST(SurfaceEditingTest, MatchSurfaceExplicitRows) { + auto plane_a = make_test_plane(); + auto plane_b = make_test_plane(); + + MatchSurfaceOptions opts; + opts.continuity = ContinuityLevel::G2; + opts.rows_to_adjust = 1; // override, only 1 row + opts.boundary_samples = 10; + + auto result = match_surface(plane_b, plane_a, BoundaryEdge::VMin, opts); + EXPECT_EQ(result.rows_adjusted, 1); +} + +// ═══════════════════════════════════════════════════════════ +// control_point_edit 测试 +// ═══════════════════════════════════════════════════════════ + +TEST(SurfaceEditingTest, ControlPointEditSingle) { + auto surf = make_test_plane(); + + std::vector> indices = {{1, 1}}; // center control point + std::vector new_pos = {{1, 1, 5}}; // lift up in z + + auto result = control_point_edit(surf, indices, new_pos); + + EXPECT_EQ(result.points_modified, 1); + EXPECT_GT(result.max_displacement, 4.0); // moved ~5 units in z + EXPECT_TRUE(result.weights_preserved); + + // Verify the control point was moved + auto modified_cp = result.modified_surface.control_points(); + EXPECT_NEAR(modified_cp[1][1].z(), 5.0, 1e-6); + // Other control points unchanged + EXPECT_NEAR(modified_cp[0][0].z(), 0.0, 1e-6); +} + +TEST(SurfaceEditingTest, ControlPointEditMultiple) { + auto surf = make_test_plane(); + + std::vector> indices = {{0, 0}, {2, 2}}; + std::vector new_pos = {{0, 0, 3}, {3, 3, 3}}; + + auto result = control_point_edit(surf, indices, new_pos); + + EXPECT_EQ(result.points_modified, 2); + EXPECT_GT(result.max_displacement, 2.0); + + auto cp = result.modified_surface.control_points(); + EXPECT_NEAR(cp[0][0].z(), 3.0, 1e-6); + EXPECT_NEAR(cp[2][2].z(), 3.0, 1e-6); + // Center unchanged + EXPECT_NEAR(cp[1][1].z(), 0.0, 1e-6); +} + +TEST(SurfaceEditingTest, ControlPointEditEmpty) { + auto surf = make_test_plane(); + + std::vector> indices; + std::vector new_pos; + + auto result = control_point_edit(surf, indices, new_pos); + + EXPECT_EQ(result.points_modified, 0); + EXPECT_NEAR(result.max_displacement, 0.0, 1e-12); + + // Should return a valid copy of original + auto cp = result.modified_surface.control_points(); + EXPECT_EQ(cp.size(), 3u); + EXPECT_NEAR(cp[0][0].z(), 0.0, 1e-6); +} + +TEST(SurfaceEditingTest, ControlPointEditPreservesKnots) { + auto surf = make_test_plane(); + + std::vector> indices = {{0, 0}}; + std::vector new_pos = {{-1, -1, 0}}; + + auto result = control_point_edit(surf, indices, new_pos); + + // Knot vectors preserved + const auto& ku = result.modified_surface.knots_u(); + const auto& kv = result.modified_surface.knots_v(); + EXPECT_EQ(ku.size(), 6u); + EXPECT_EQ(kv.size(), 6u); + EXPECT_EQ(result.modified_surface.degree_u(), 2); + EXPECT_EQ(result.modified_surface.degree_v(), 2); +} + +} // namespace test +} // namespace vde::curves