From 73df04d5cb722bb6d7b18378b511cc9773f4df5b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E8=8C=82=E4=B9=8B=E9=92=B3?= Date: Sun, 26 Jul 2026 20:58:31 +0800 Subject: [PATCH] feat(v5-M2): G2/G3 continuity + surface analysis + extension + N-side fill + advanced blend MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit M2.1 — G2/G3 连续性分析 (Agent #0): - surface_continuity.h/.cpp: G0/G1/G2/G3 curve/surface detection - Weingarten equation for curvature, Frénet frame, zebra stripe - surface_analysis.h/.cpp: curvature_map, deviation_analysis, curvature_comb - 24 tests (12 continuity + 12 analysis) M2.2 — 曲面延伸 + N边填充 (Agent #1): - surface_extension.h/.cpp: extend_surface(G1/G2), n_sided_fill, blend_surfaces - Coons patch generalization for N-sided holes - 16/16 tests passed in Docker container M2.3 — 高级过渡曲面 (Agent #2): - advanced_blend.h/.cpp: real implementations replacing stubs - variable_radius_blend, multi_face_blend, rolling_ball_blend, face_face_blend - Ball-rolling envelope + corner sphere filling - 15+ tests with validate() verification --- include/vde/brep/advanced_blend.h | 126 +++- include/vde/curves/surface_analysis.h | 228 ++++++ include/vde/curves/surface_continuity.h | 155 ++++ include/vde/curves/surface_extension.h | 116 +++ src/CMakeLists.txt | 3 + src/brep/advanced_blend.cpp | 917 ++++++++++++++++++++++- src/curves/surface_analysis.cpp | 537 +++++++++++++ src/curves/surface_continuity.cpp | 610 +++++++++++++++ src/curves/surface_extension.cpp | 636 ++++++++++++++++ tests/brep/CMakeLists.txt | 1 + tests/brep/test_advanced_blend.cpp | 158 ++++ tests/curves/CMakeLists.txt | 3 + tests/curves/test_surface_analysis.cpp | 181 +++++ tests/curves/test_surface_continuity.cpp | 195 +++++ tests/curves/test_surface_extension.cpp | 303 ++++++++ 15 files changed, 4158 insertions(+), 11 deletions(-) create mode 100644 include/vde/curves/surface_analysis.h create mode 100644 include/vde/curves/surface_continuity.h create mode 100644 include/vde/curves/surface_extension.h create mode 100644 src/curves/surface_analysis.cpp create mode 100644 src/curves/surface_continuity.cpp create mode 100644 src/curves/surface_extension.cpp create mode 100644 tests/brep/test_advanced_blend.cpp create mode 100644 tests/curves/test_surface_analysis.cpp create mode 100644 tests/curves/test_surface_continuity.cpp create mode 100644 tests/curves/test_surface_extension.cpp diff --git a/include/vde/brep/advanced_blend.h b/include/vde/brep/advanced_blend.h index 4525b4d..4c7dcbd 100644 --- a/include/vde/brep/advanced_blend.h +++ b/include/vde/brep/advanced_blend.h @@ -1,8 +1,126 @@ #pragma once +/** + * @file advanced_blend.h + * @brief 高级过渡曲面（Advanced Blend）操作 + * + * 提供多种超出简单恒定半径圆角的过渡曲面算法： + * - 沿边变半径过渡（variable_radius_blend） + * - 多面同时过渡（multi_face_blend） + * - 滚动球算法（rolling_ball_blend） + * - 两面间恒定半径过渡（face_face_blend） + * + * ## 与 fillet/fillet_variable 的关系 + * + * modeling.h 中的 fillet 和 fillet_variable 是面向最终用户的高层 API， + * 本模块提供更底层的过渡曲面算法实现。两者可配合使用。 + * + * @ingroup brep + */ #include "vde/brep/brep.h" #include + namespace vde::brep { -[[nodiscard]] BrepModel variable_radius_blend(const BrepModel& body, int edge_id, double r_start, double r_end); -[[nodiscard]] BrepModel multi_face_blend(const BrepModel& body, const std::vector& face_ids, double radius); -[[nodiscard]] BrepModel rolling_ball_blend(const BrepModel& body, int edge_id, double radius); -} // namespace vde::brep + +/** + * @brief 沿边变半径过渡 + * + * 沿指定边的参数方向，在 r_start 和 r_end 之间线性插值半径， + * 每段独立构建过渡曲面片段，最终拼接为完整过渡。 + * + * 算法： + * 1. 定位目标边及其两个邻面 + * 2. 计算两面法向量及角平分线方向 + * 3. 沿边曲线采样 N 个截面，每个截面处计算半径和切点 + * 4. 每相邻两截面间构建一张过渡曲面（NURBS quad patch） + * 5. 重建两个邻面，切除过渡区域 + * 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 samples 越大截面过渡越光滑，但面和边数增加 + */ +[[nodiscard]] BrepModel variable_radius_blend( + const BrepModel& body, int edge_id, + double r_start, double r_end, int samples = 16); + +/** + * @brief 多面同时过渡 + * + * 对多个面同时施加恒定半径的过渡处理。每个目标面独立向内偏移， + * 目标面之间的交角处构造球面填充，非目标面保持不变。 + * + * 算法： + * 1. 对每个目标面，沿面法向量向内偏移半径距离 + * 2. 检测目标面之间的共享边，在交线处构造圆柱过渡 + * 3. 检测三面以上交角顶点，用球面贴片填充 + * 4. 复制非目标面，组装壳和体 + * + * @param body 输入实体 + * @param face_ids 目标面索引列表 + * @param radius 过渡半径（> 0） + * @return 过渡后的新实体 + * + * @pre face_ids 非空，所有索引有效，radius > 0 + */ +[[nodiscard]] BrepModel multi_face_blend( + const BrepModel& body, const std::vector& face_ids, double radius); + +/** + * @brief 滚动球算法过渡 + * + * 模拟一个半径为 radius 的球体沿目标边滚动，球体的扫掠包络面即为过渡曲面。 + * 球体始终保持与边的两个邻面接触，其中心沿边的等距曲线移动。 + * + * 算法： + * 1. 计算滚动球中心轨迹（沿边的 offset 曲线，offset = r / cos(half_angle)） + * 2. 沿边采样球心位置 + * 3. 每个球心位置计算球与两面的切点 + * 4. 相邻截面间构建过渡曲面 patch（近似球体扫掠包络） + * 5. 重建邻面、复制非影响面、组装壳和体 + * + * 对直边退化为圆柱过渡，对曲边产生自然的变曲率过渡。 + * + * @param body 输入实体 + * @param edge_id 目标边索引 + * @param radius 球体半径（> 0） + * @param samples 沿边的球心采样数（默认 20） + * @return 过渡后的新实体 + * + * @pre edge_id 有效，radius > 0 + */ +[[nodiscard]] BrepModel rolling_ball_blend( + const BrepModel& body, int edge_id, + double radius, int samples = 20); + +/** + * @brief 两面间恒定半径过渡 + * + * 在两个任意面之间构造恒定半径的过渡曲面。与单边过渡不同， + * 此函数处理两个面之间的整个共享或邻近区域。 + * + * 算法： + * 1. 找到两面之间的共享边，或确定两面的邻近区域 + * 2. 沿共享边界计算过渡曲面（恒定半径圆柱面） + * 3. 切掉两面在过渡区域的部分 + * 4. 插入过渡曲面，组装壳和体 + * + * @param body 输入实体 + * @param face_a 第一个面索引 + * @param face_b 第二个面索引 + * @param radius 过渡半径（> 0） + * @param samples 沿边界的截面采样数（默认 16） + * @return 过渡后的新实体 + * + * @pre face_a 和 face_b 有效，radius > 0 + */ +[[nodiscard]] BrepModel face_face_blend( + const BrepModel& body, int face_a, int face_b, + double radius, int samples = 16); + +} // namespace vde::brep diff --git a/include/vde/curves/surface_analysis.h b/include/vde/curves/surface_analysis.h new file mode 100644 index 0000000..04d37c2 --- /dev/null +++ b/include/vde/curves/surface_analysis.h @@ -0,0 +1,228 @@ +#pragma once +/** + * @file surface_analysis.h + * @brief 曲面分析工具集 — 曲率映射、斑马纹、曲率梳、偏差分析、拔模角 + * + * 提供 CAD 级曲面质量分析功能： + * - 高斯/平均/主曲率网格 + * - 斑马纹反射模拟 + * - 曲率梳可视化数据 + * - 曲面间偏差分析 + * - 拔模角计算（增强现有实现） + * + * @ingroup curves + */ + +#include "vde/curves/nurbs_curve.h" +#include "vde/curves/nurbs_surface.h" +#include "vde/core/point.h" +#include "vde/brep/brep.h" +#include +#include +#include +#include + +namespace vde::curves { + +using core::Point3D; +using core::Vector3D; + +// ═══════════════════════════════════════════════════════════ +// Curvature Map Types +// ═══════════════════════════════════════════════════════════ + +/// 单点的曲率数据 +struct CurvaturePoint { + double gaussian = 0.0; ///< 高斯曲率 K = κ₁·κ₂ + double mean = 0.0; ///< 平均曲率 H = (κ₁+κ₂)/2 + double k1 = 0.0; ///< 最大主曲率 κ₁ + double k2 = 0.0; ///< 最小主曲率 κ₂ + Vector3D d1; ///< 最大主方向 + Vector3D d2; ///< 最小主方向 +}; + +/// 曲率网格（res_u+1 × res_v+1） +struct CurvatureMap { + int res_u, res_v; ///< 网格分辨率 + std::vector> grid; ///< grid[i][j] 对应 (u_i, v_j) + + /// 高斯曲率极值 + double min_gaussian = std::numeric_limits::max(); + double max_gaussian = -std::numeric_limits::max(); + + /// 平均曲率极值 + double min_mean = std::numeric_limits::max(); + double max_mean = -std::numeric_limits::max(); + + /// 是否为可展曲面（高斯曲率 ≈ 0） + [[nodiscard]] bool is_developable(double tol = 1e-9) const; +}; + +/** + * @brief 计算 NURBS 曲面的曲率映射 + * + * 在参数域均匀采样，计算每个采样点的高斯曲率、平均曲率、 + * 主曲率及主方向。 + * + * @param surf NURBS 曲面 + * @param res_u u 方向分辨率（单元格数） + * @param res_v v 方向分辨率（单元格数） + * @return CurvatureMap + */ +[[nodiscard]] CurvatureMap curvature_map( + const NurbsSurface& surf, + int res_u, int res_v); + +// ═══════════════════════════════════════════════════════════ +// Zebra Stripe +// ═══════════════════════════════════════════════════════════ + +/// 斑马纹反射数据 +struct ZebraStripeResult { + int res_u, res_v; ///< 网格分辨率 + std::vector> intensity; ///< intensity[i][j] ∈ [0,1] + Vector3D light_direction; ///< 光照方向（归一化） +}; + +/** + * @brief 计算斑马纹反射数据（等照度线模拟） + * + * 模拟等间距平行光源在曲面上的反射效果，用于检测曲面光顺性。 + * 斑马纹的密度和均匀度反映曲面的曲率变化。 + * + * @param surf NURBS 曲面 + * @param light_dir 光照方向（自动归一化） + * @param res_u u 方向分辨率 + * @param res_v v 方向分辨率 + * @return 反射强度网格 + */ +[[nodiscard]] ZebraStripeResult zebra_stripe( + const NurbsSurface& surf, + const Vector3D& light_dir, + int res_u, int res_v); + +// ═══════════════════════════════════════════════════════════ +// Curvature Comb +// ═══════════════════════════════════════════════════════════ + +/// 曲率梳数据点 +struct CurvatureCombPoint { + Point3D curve_point; ///< 曲线上的点 + Vector3D curvature_vec; ///< 曲率向量（方向=主法线，长度∝曲率） + double curvature; ///< 曲率大小 +}; + +/// 曲率梳完整数据 +struct CurvatureComb { + std::vector points; + double min_curvature = std::numeric_limits::max(); + double max_curvature = -std::numeric_limits::max(); + double scale = 1.0; ///< 缩放因子 +}; + +/** + * @brief 计算曲线曲率梳数据 + * + * 在曲线上均匀采样，计算各点的曲率向量。 + * 曲率梳将曲率大小可视化为从曲线向外延伸的线段长度。 + * + * @param curve NURBS 曲线 + * @param res 采样分辨率（点数） + * @param scale 曲率缩放因子（默认自动适配） + * @return 曲率梳数据 + */ +[[nodiscard]] CurvatureComb curvature_comb( + const NurbsCurve& curve, + int res, + double scale = 0.0); + +// ═══════════════════════════════════════════════════════════ +// Deviation Analysis +// ═══════════════════════════════════════════════════════════ + +/// 偏差分析结果 +struct DeviationResult { + double max_positive = 0.0; ///< 最大正偏差（surf_b 在 surf_a 外侧） + double max_negative = 0.0; ///< 最大负偏差（surf_b 在 surf_a 内侧） + double rms = 0.0; ///< 均方根偏差 + double mean_absolute = 0.0; ///< 平均绝对偏差 + size_t sample_count = 0; ///< 采样点数量 + std::vector deviations; ///< 各采样点的有符号偏差 + std::vector sample_points_a; ///< surf_a 上的采样点 + std::vector sample_points_b; ///< surf_b 上对应的最近点 +}; + +/** + * @brief 计算两曲面的偏差分析 + * + * 在 surf_a 参数域均匀采样，对每个采样点寻找 surf_b 上的最近点， + * 沿 surf_a 法线方向计算有符号偏差。 + * + * @param surf_a 基准曲面 + * @param surf_b 对比曲面 + * @param samples 采样点数量（总点数 = samples×samples） + * @return 偏差分析结果 + */ +[[nodiscard]] DeviationResult deviation_analysis( + const NurbsSurface& surf_a, + const NurbsSurface& surf_b, + int samples); + +// ═══════════════════════════════════════════════════════════ +// Draft Face Angle (Enhanced) +// ═══════════════════════════════════════════════════════════ + +/** + * @brief 计算 NURBS 曲面上指定参数点的拔模角 + * + * 拔模角 = angle between surface normal and pull direction minus 90° + * 正值表示面沿拔模方向向外张开。 + * + * @param surf NURBS 曲面 + * @param u u 参数 + * @param v v 参数 + * @param pull_dir 拔模方向（自动归一化） + * @return 拔模角（弧度） + */ +[[nodiscard]] double draft_face_angle( + const NurbsSurface& surf, + double u, double v, + const Vector3D& pull_dir); + +/** + * @brief 计算 B-Rep 面的拔模角（增强版） + * + * 在面的参数域中采样多个点，返回面积加权平均拔模角。 + * + * @param body B-Rep 模型 + * @param face_id 面 ID + * @param pull_dir 拔模方向 + * @param samples 每方向采样点数（默认 10） + * @return 面积加权平均拔模角（弧度），若面无效则返回 0 + */ +[[nodiscard]] double draft_face_angle( + const brep::BrepModel& body, + int face_id, + const Vector3D& pull_dir, + int samples = 10); + +/** + * @brief 计算 B-Rep 面的拔模角分布 + * + * 返回面上多个采样点的拔模角数据。 + */ +struct DraftAngleDistribution { + double min_angle, max_angle; ///< 拔模角范围 + double mean_angle; ///< 平均拔模角 + double area_weighted_angle; ///< 面积加权平均 + std::vector sample_angles; ///< 各采样点拔模角 + std::vector> sample_params; ///< 对应参数 (u,v) +}; + +[[nodiscard]] DraftAngleDistribution draft_face_angle_distribution( + const brep::BrepModel& body, + int face_id, + const Vector3D& pull_dir, + int samples = 15); + +} // namespace vde::curves diff --git a/include/vde/curves/surface_continuity.h b/include/vde/curves/surface_continuity.h new file mode 100644 index 0000000..1b87bec --- /dev/null +++ b/include/vde/curves/surface_continuity.h @@ -0,0 +1,155 @@ +#pragma once +/** + * @file surface_continuity.h + * @brief 曲面连续性分析 — G0/G1/G2/G3 连续性检测 + * + * 提供曲线间和曲面间沿公共边的连续性分类。 + * + * ### 连续性等级 + * + * | 等级 | 条件 | + * |------|------| + * | G0 | 位置连续：两端点重合 | + * | G1 | 切向连续：一阶导数共线（方向一致） | + * | G2 | 曲率连续：曲率大小和方向匹配 | + * | G3 | 曲率变化率连续：曲率的变化率匹配 | + * + * @ingroup curves + */ + +#include "vde/curves/nurbs_curve.h" +#include "vde/curves/nurbs_surface.h" +#include "vde/core/point.h" +#include +#include + +namespace vde::curves { + +// ═══════════════════════════════════════════════════════════ +// Continuity Level Enum +// ═══════════════════════════════════════════════════════════ + +/// 连续性等级 +enum class ContinuityLevel { + None, ///< 不连续（甚至不满足 G0） + G0, ///< 位置连续 + G1, ///< 切向连续 + G2, ///< 曲率连续 + G3, ///< 曲率变化率连续 +}; + +/// 将 ContinuityLevel 转为可读字符串 +[[nodiscard]] const char* to_string(ContinuityLevel level); + +// ═══════════════════════════════════════════════════════════ +// Curve Continuity +// ═══════════════════════════════════════════════════════════ + +/** + * @brief 分析两条 NURBS 曲线在公共参数点的连续性 + * + * 在给定参数点计算两曲线的位置、一阶、二阶、三阶导数， + * 按 G0→G1→G2→G3 逐级判定。 + * + * @param curve_a 第一条 NURBS 曲线 + * @param t_a curve_a 上公共点的参数值 + * @param curve_b 第二条 NURBS 曲线 + * @param t_b curve_b 上公共点的参数值 + * @param tol 位置容差（默认 1e-6） + * @return 满足的最高连续性等级 + * + * @note G1 要求一阶导数共线（方向相同），G2 要求曲率向量一致， + * G3 要求曲率变化率一致。判定是累积的：G3 自动满足 G2/G1/G0。 + */ +[[nodiscard]] ContinuityLevel continuity_type( + const NurbsCurve& curve_a, double t_a, + const NurbsCurve& curve_b, double t_b, + double tol = 1e-6); + +/** + * @brief 重载：在两条曲线端点处检测连续性 + * + * 自动探测 curve_a 和 curve_b 各两个端点，检查是否有重合端点对， + * 然后调用上述函数分析该连接点的连续性等级。 + * + * @param curve_a 第一条曲线 + * @param curve_b 第二条曲线 + * @param tol 位置容差 + * @return 若存在重合端点则返回连续性等级，否则返回 ContinuityLevel::None + */ +[[nodiscard]] ContinuityLevel continuity_type( + const NurbsCurve& curve_a, + const NurbsCurve& curve_b, + double tol = 1e-6); + +// ═══════════════════════════════════════════════════════════ +// Surface Continuity +// ═══════════════════════════════════════════════════════════ + +/** + * @brief 沿公共等参边分析两张 NURBS 曲面的连续性 + * + * 在沿公共边的采样点处评估曲面连续性，返回最弱点等级。 + * 例如 surf_a 的 u=u_max 边与 surf_b 的 u=u_min 边相邻。 + * + * @param surf_a 第一张曲面 + * @param surf_b 第二张曲面 + * @param edge_param_range_a 公共边在 surf_a 上的参数范围 {{u, v_start}, {u, v_end}} + * @param edge_param_range_b 公共边在 surf_b 上的参数范围 {{u, v_start}, {u, v_end}} + * @param samples 沿边的采样点数量（默认 20） + * @param tol 位置容差（默认 1e-6） + * @return 所有采样点中的最低连续性等级 + * + * @note 若曲面在采样点间连续但曲面方向（法线）相反， + * 仍判定为 G0（因为 G1 要求一阶导数方向一致）。 + */ +[[nodiscard]] ContinuityLevel surface_continuity( + const NurbsSurface& surf_a, + const NurbsSurface& surf_b, + std::pair, std::pair> edge_param_range_a, + std::pair, std::pair> edge_param_range_b, + int samples = 20, + double tol = 1e-6); + +/** + * @brief 沿公共边自动检测并分析两张曲面连续性 + * + * 自动检测 surf_a 和 surf_b 的公共边（通过比较四条等参边界）， + * 然后调用 surface_continuity 分析。 + * + * @param surf_a 第一张曲面 + * @param surf_b 第二张曲面 + * @param samples 采样点数量 + * @param tol 位置容差 + * @return 所有公共边中的最低连续性等级（若无公共边则返回 None） + */ +[[nodiscard]] ContinuityLevel surface_continuity( + const NurbsSurface& surf_a, + const NurbsSurface& surf_b, + int samples = 20, + double tol = 1e-6); + +/** + * @brief 曲面间连续性详细报告 + * + * 返回每个采样点的连续性等级向量，便于诊断。 + */ +struct SurfaceContinuityReport { + ContinuityLevel worst; ///< 最弱等级 + ContinuityLevel best; ///< 最佳等级 + double g1_ratio; ///< G1 通过率 (0~1) + double g2_ratio; ///< G2 通过率 (0~1) + std::vector per_sample; ///< 各采样点等级 + std::vector position_errors; ///< 各采样点位置误差 +}; + +/** + * @brief 获取曲面连续性详细报告 + */ +[[nodiscard]] SurfaceContinuityReport surface_continuity_report( + const NurbsSurface& surf_a, + const NurbsSurface& surf_b, + int samples = 50, + double tol = 1e-6); + +} // namespace vde::curves diff --git a/include/vde/curves/surface_extension.h b/include/vde/curves/surface_extension.h new file mode 100644 index 0000000..c182041 --- /dev/null +++ b/include/vde/curves/surface_extension.h @@ -0,0 +1,116 @@ +#pragma once +#include "vde/curves/nurbs_surface.h" +#include "vde/curves/nurbs_curve.h" +#include "vde/core/point.h" +#include +#include + +namespace vde::curves { + +using core::Point3D; +using core::Vector3D; + +/** + * @brief 连续性级别 + * @ingroup curves + */ +enum class Continuity { + C0 = 0, ///< 仅位置连续 + G1 = 1, ///< 切线/切平面连续 + G2 = 2, ///< 曲率连续 +}; + +/** + * @brief 曲面延伸 + * + * 沿等参边界延伸曲面，支持 G1 切平面连续。 + * 对于 G2 连续，在约束方向构造 3 排控制点以保持曲率。 + * + * @param surface 输入 NURBS 曲面 + * @param edge_id 边界边索引（0=umin, 1=umax, 2=vmin, 3=vmax） + * @param distance 延伸距离（参数空间） + * @param continuity 期望的连续性级别（默认 G1） + * @return 延伸后的 NURBS 曲面 + * + * @note G1 模式：追加 1 排控制点，沿切向延伸 + * @note G2 模式：追加 2 排控制点，第 2 排由曲率约束计算 + * @note 延伸仅沿切向线性/二次外推，不保证与原曲面精确 G2 连续 + * + * @code{.cpp} + * auto ext = extend_surface(surf, 1, 0.5, Continuity::G1); + * @endcode + * @ingroup curves + */ +[[nodiscard]] NurbsSurface extend_surface(const NurbsSurface& surface, + int edge_id, double distance, + Continuity continuity = Continuity::G1); + +/** + * @brief N 边孔洞填充 + * + * 用 N 条边界曲线构造 Coons patch 泛化的填充曲面。 + * 支持三角形(N=3)、四边形(N=4)、五边形(N=5)等任意边数。 + * + * 算法：将边界曲线映射到单位圆参数域，使用径向基函数插值 + * 构造中心点，再以双线性 Coons 方式生成内部控制点。 + * + * @param boundary_curves N 条边界曲线，需围成闭合环（首尾相连） + * @param continuity 期望的连续性级别（默认 G2） + * @return 填充 NURBS 曲面 + * + * @note G1 模式：单层中间环，线性插值 + * @note G2 模式：双层中间环，保持曲率衰减 + * @note 至少需要 3 条边界曲线 + * + * @code{.cpp} + * auto fill = n_sided_fill({c1, c2, c3}, Continuity::G2); // 三角形孔 + * @endcode + * @ingroup curves + */ +[[nodiscard]] NurbsSurface n_sided_fill(const std::vector& boundary_curves, + Continuity continuity = Continuity::G2); + +/** + * @brief 两面间过渡曲面 + * + * 检测两面公共边界，沿该边构造圆弧或 S 形过渡曲面。 + * 与 nurbs_operations 中的 blend_surfaces 不同，此版本自动检测公共边。 + * + * 算法： + * 1. 检测 surf_a 和 surf_b 的共享边界 + * 2. 沿共享边界两侧各偏移 radius 提取曲线 + * 3. 在两条偏移曲线间构造过渡曲面（直纹面或 Coons 面） + * + * @param surf_a 第一个曲面 + * @param surf_b 第二个曲面 + * @param radius 过渡半径 + * @return 过渡 NURBS 曲面 + * + * @note 如果两面无公共边界，返回退化的空曲面 + * @note 过渡面在边界处与原曲面 G1 连续（等参方向） + * + * @code{.cpp} + * auto blend = blend_surfaces(surf_a, surf_b, 0.3); + * @endcode + * @ingroup curves + */ +[[nodiscard]] NurbsSurface blend_surfaces(const NurbsSurface& surf_a, + const NurbsSurface& surf_b, + double radius); + +/** + * @brief 检查两曲面沿共享边界的 G1 连续性 + * + * @param surf_a 第一个曲面 + * @param surf_b 第二个曲面 + * @param edge surf_a 的边界边索引（0=umin, 1=umax, 2=vmin, 3=vmax） + * @param samples 沿边界的采样点数（默认 10） + * @return 所有采样点均满足 G1 连续时返回 true + * + * @ingroup curves + */ +[[nodiscard]] bool is_g1_continuous(const NurbsSurface& surf_a, + const NurbsSurface& surf_b, + int edge, int samples = 10); + +} // namespace vde::curves diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 85c0f7e..d0a074f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -52,6 +52,9 @@ add_library(vde_curves STATIC curves/parallel_intersection.cpp curves/surface_fairing.cpp curves/exact_offset.cpp + curves/surface_extension.cpp + curves/surface_continuity.cpp + curves/surface_analysis.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 9dc56fb..11cd71d 100644 --- a/src/brep/advanced_blend.cpp +++ b/src/brep/advanced_blend.cpp @@ -1,12 +1,915 @@ #include "vde/brep/advanced_blend.h" +#include "vde/brep/modeling.h" +#include "vde/brep/brep_validate.h" +#include "vde/curves/nurbs_curve.h" +#include "vde/curves/nurbs_surface.h" +#include +#include +#include +#include +#include +#include + namespace vde::brep { -BrepModel variable_radius_blend(const BrepModel& body, int edge_id, double r_start, double r_end) { - (void)edge_id;(void)r_start;(void)r_end; return body; +namespace { + +// ═══════════════════════════════════════════════ +// Helper: geometry utilities +// ═══════════════════════════════════════════════ + +using core::Point3D; +using core::Vector3D; +using core::AABB3D; +using curves::NurbsCurve; +using curves::NurbsSurface; + +/// Compute face normal from its first 3 vertices +Vector3D face_normal(const BrepModel& body, int face_idx) { + auto es = body.face_edges(face_idx); + if (es.size() < 3) return Vector3D::UnitZ(); + const auto& e0 = body.edge(es[0]); + const auto& e1 = body.edge(es[1]); + Point3D p0 = body.vertex_by_id(e0.v_start).point; + Point3D p1 = body.vertex_by_id(e0.v_end).point; + Point3D p2 = body.vertex_by_id(e1.v_end).point; + Vector3D n = (p1 - p0).cross(p2 - p0); + double len = n.norm(); + if (len > 1e-12) return n / len; + return Vector3D::UnitZ(); } -BrepModel multi_face_blend(const BrepModel& body, const std::vector& face_ids, double radius) { - (void)face_ids;(void)radius; return body; + +/// Linear interpolation of point +Point3D lerp(const Point3D& a, const Point3D& b, double t) { + return { + a.x() + (b.x() - a.x()) * t, + a.y() + (b.y() - a.y()) * t, + a.z() + (b.z() - a.z()) * t + }; } -BrepModel rolling_ball_blend(const BrepModel& body, int edge_id, double radius) { - (void)edge_id;(void)radius; return body; + +/// Build a degree-1×1 plane surface from 4 corner points +NurbsSurface make_plane(const Point3D& p00, const Point3D& p10, + const Point3D& p11, const Point3D& p01) { + std::vector> grid = {{p00, p01}, {p10, p11}}; + return NurbsSurface(grid, {0, 0, 1, 1}, {0, 0, 1, 1}, {}, 1, 1); } -} // namespace vde::brep + +/// Build degree-2×2 arc patch for blend surface segment +NurbsSurface make_blend_patch( + const Point3D& t1a, const Point3D& ma, const Point3D& t2a, + const Point3D& t1b, const Point3D& mb, const Point3D& t2b) { + std::vector> grid = { + {t1a, ma, t2a}, + {t1b, mb, t2b} + }; + return NurbsSurface(grid, + {0, 0, 1, 1}, + {0, 0, 0, 1, 1, 1}, + {}, 1, 2); +} + +/// Build a degree-2×2 spherical corner patch +NurbsSurface make_sphere_corner_patch( + const Point3D& corner, const Vector3D& n1, const Vector3D& n2, + const Vector3D& n3, double r) { + // 3×3 control grid approximating spherical octant at corner + std::vector> grid(3); + for (int i = 0; i < 3; ++i) { + double u = i / 2.0; + for (int j = 0; j < 3; ++j) { + double v = j / 2.0; + // Spherical interpolation between the 3 tangent directions + Vector3D dir = (n1 * (1 - u) + n2 * u) * (1 - v) + n3 * v; + double dlen = dir.norm(); + if (dlen < 1e-12) dir = n1; + else dir = dir / dlen; + grid[i].push_back({ + corner.x() + dir.x() * r, + corner.y() + dir.y() * r, + corner.z() + dir.z() * r + }); + } + } + std::vector ku = {0, 0, 0, 1, 1, 1}; + std::vector kv = {0, 0, 0, 1, 1, 1}; + std::vector> w(3, std::vector(3, 1.0)); + return NurbsSurface(grid, ku, kv, w, 2, 2); +} + +/// Find face2 that shares the given edge with face1 +int find_other_face(const BrepModel& body, int edge_id, int face1) { + auto efs = body.edge_faces(edge_id); + for (int fi : efs) { + if (fi != face1) return fi; + } + return -1; +} + +/// Check if two edges are the same (matching endpoints in either order) +bool edges_same(const BrepModel& body, int e1, int e2) { + const auto& a = body.edge(e1); + const auto& b = body.edge(e2); + Point3D a0 = body.vertex_by_id(a.v_start).point; + Point3D a1 = body.vertex_by_id(a.v_end).point; + Point3D b0 = body.vertex_by_id(b.v_start).point; + Point3D b1 = body.vertex_by_id(b.v_end).point; + double eps = 1e-7; + return ((a0 - b0).norm() < eps && (a1 - b1).norm() < eps) || + ((a0 - b1).norm() < eps && (a1 - b0).norm() < eps); +} + +/// Get the curve for an edge, or construct a line if no curve +NurbsCurve edge_curve(const BrepModel& body, int edge_id) { + const auto& e = body.edge(edge_id); + if (e.curve) return *e.curve; + Point3D p0 = body.vertex_by_id(e.v_start).point; + Point3D p1 = body.vertex_by_id(e.v_end).point; + return NurbsCurve({p0, p1}, {0, 0, 1, 1}, {1, 1}, 1); +} + +// ═══════════════════════════════════════════════ +// Helper: build quad face in result model +// ═══════════════════════════════════════════════ +int add_quad(BrepModel& m, const Point3D& v0, const Point3D& v1, + const Point3D& v2, const Point3D& v3, const NurbsSurface& surf) { + int iv0 = m.add_vertex(v0), iv1 = m.add_vertex(v1); + int iv2 = m.add_vertex(v2), iv3 = m.add_vertex(v3); + int e01 = m.add_edge(iv0, iv1), e12 = m.add_edge(iv1, iv2); + int e23 = m.add_edge(iv2, iv3), e30 = m.add_edge(iv3, iv0); + int sid = m.add_surface(surf); + return m.add_face(sid, {m.add_loop({e01, e12, e23, e30}, true)}); +} + +/// Add a triangular face +int add_tri(BrepModel& m, const Point3D& v0, const Point3D& v1, + const Point3D& v2, const NurbsSurface& surf) { + int iv0 = m.add_vertex(v0), iv1 = m.add_vertex(v1), iv2 = m.add_vertex(v2); + int e01 = m.add_edge(iv0, iv1), e12 = m.add_edge(iv1, iv2), + e20 = m.add_edge(iv2, iv0); + int sid = m.add_surface(surf); + return m.add_face(sid, {m.add_loop({e01, e12, e20}, true)}); +} + +/// Copy an unaffected face from source to result model +void copy_face(const BrepModel& src, int face_idx, BrepModel& dst, + std::map& vertex_remap) { + const auto& f_src = src.face(face_idx); + for (int li : f_src.loops) { + const auto& lop = src.loop_by_id(li); + std::vector new_es; + for (int ei : lop.edges) { + const auto& e_src = src.edge(ei); + int sv = e_src.v_start, ev = e_src.v_end; + if (!vertex_remap.count(sv)) + vertex_remap[sv] = dst.add_vertex(src.vertex_by_id(sv).point); + if (!vertex_remap.count(ev)) + vertex_remap[ev] = dst.add_vertex(src.vertex_by_id(ev).point); + new_es.push_back(dst.add_edge(vertex_remap[sv], vertex_remap[ev])); + } + int new_lp = dst.add_loop(new_es, lop.is_outer); + int new_surf = dst.add_surface(src.surface(f_src.surface_id)); + dst.add_face(new_surf, {new_lp}); + } +} + +/// Finish building the result model: wrap faces in a shell and body +void finish_model(BrepModel& result, const std::string& name) { + std::vector all_faces; + for (size_t fi = 0; fi < result.num_faces(); ++fi) + all_faces.push_back(result.face(static_cast(fi)).id); + if (!all_faces.empty()) { + int sh = result.add_shell(all_faces, true); + result.add_body({sh}, name); + } +} + +} // anonymous namespace + +// ═══════════════════════════════════════════════ +// variable_radius_blend +// ═══════════════════════════════════════════════ + +BrepModel variable_radius_blend(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 + 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() + }; + + // Tangent points on each face + 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; + } + + // Edge direction for perpendicular computation + Vector3D edge_dir = (edge_pts.back().x() - edge_pts[0].x() > 0 || + edge_pts.back().y() - edge_pts[0].y() > 0 || + edge_pts.back().z() - edge_pts[0].z() > 0) ? + Vector3D{edge_pts.back().x() - edge_pts[0].x(), + edge_pts.back().y() - edge_pts[0].y(), + edge_pts.back().z() - edge_pts[0].z()} : Vector3D::UnitX(); + double ed_len = edge_dir.norm(); + if (ed_len > 1e-12) edge_dir = edge_dir / ed_len; + + BrepModel result; + + // ── Rebuild face_a: follow loop order, replace blend 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); + // Find blend edge position in the loop + 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; + + // Build new edges in cyclic order, starting right after blend edge + std::vector new_es; + // Insert the new blend tangent edge + int nv0 = result.add_vertex(p_first); + int nv1 = result.add_vertex(p_last); + new_es.push_back(result.add_edge(nv0, nv1)); + + // Follow remaining edges in loop order + 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; // only process the first loop (outer) + } + } + + // ── Rebuild face_b: follow loop order (reversed orientation) ── + { + 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; + // Reversed blend edge for opposite face + 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 ── + 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); + } + + finish_model(result, "variable_radius_blend"); + return result; +} + +// ═══════════════════════════════════════════════ +// multi_face_blend +// ═══════════════════════════════════════════════ + +BrepModel multi_face_blend(const BrepModel& body, + const std::vector& face_ids, double radius) { + if (face_ids.empty() || radius <= 0.0) return body; + + std::set target_set(face_ids.begin(), face_ids.end()); + + // ── Gather all vertices of target faces ── + std::set target_vertices; // vertex IDs on target faces + std::map vertex_normal; // accumulated offset direction per vertex + std::map vertex_pos; // original position + + for (int fi : face_ids) { + auto f_edges = body.face_edges(fi); + Vector3D fn = face_normal(body, fi); + for (int ei : f_edges) { + const auto& e = body.edge(ei); + target_vertices.insert(e.v_start); + target_vertices.insert(e.v_end); + vertex_pos[e.v_start] = body.vertex_by_id(e.v_start).point; + vertex_pos[e.v_end] = body.vertex_by_id(e.v_end).point; + // Accumulate face normals for vertex offset direction + if (!vertex_normal.count(e.v_start)) + vertex_normal[e.v_start] = Vector3D{0,0,0}; + if (!vertex_normal.count(e.v_end)) + vertex_normal[e.v_end] = Vector3D{0,0,0}; + vertex_normal[e.v_start] = { + vertex_normal[e.v_start].x() + fn.x(), + vertex_normal[e.v_start].y() + fn.y(), + vertex_normal[e.v_start].z() + fn.z() + }; + vertex_normal[e.v_end] = { + vertex_normal[e.v_end].x() + fn.x(), + vertex_normal[e.v_end].y() + fn.y(), + vertex_normal[e.v_end].z() + fn.z() + }; + } + } + + // Normalize vertex offset directions + for (auto& [vid, dir] : vertex_normal) { + double len = dir.norm(); + if (len > 1e-12) dir = dir / len; + } + + BrepModel result; + + // ── Build offset vertices ── + std::map offset_vtx; // original vertex ID → new vertex index + for (int vid : target_vertices) { + auto dir = vertex_normal[vid]; + auto pos = vertex_pos[vid]; + Point3D offset_pt = { + pos.x() - dir.x() * radius, + pos.y() - dir.y() * radius, + pos.z() - dir.z() * radius + }; + offset_vtx[vid] = result.add_vertex(offset_pt); + } + + // ── Rebuild each target face with offset vertices ── + for (int fi : face_ids) { + const auto& f_src = body.face(fi); + std::vector new_es; + for (int li : f_src.loops) { + const auto& lop = body.loop_by_id(li); + for (int ei : lop.edges) { + const auto& e_src = body.edge(ei); + int nv0 = offset_vtx[e_src.v_start]; + int nv1 = offset_vtx[e_src.v_end]; + new_es.push_back(result.add_edge(nv0, nv1)); + } + } + int new_lp = result.add_loop(new_es, true); + int new_surf = result.add_surface(body.surface(f_src.surface_id)); + result.add_face(new_surf, {new_lp}); + } + + // ── Copy non-target faces (with offset vertices where they touch target faces) ── + std::map normal_vremap; // original vertex ID → result vertex + for (size_t fi = 0; fi < body.num_faces(); ++fi) { + int fidx = static_cast(fi); + if (target_set.count(fidx)) continue; + + const auto& f_src = body.face(fidx); + for (int li : f_src.loops) { + const auto& lop = body.loop_by_id(li); + std::vector new_es; + for (int ei : lop.edges) { + const auto& e_src = body.edge(ei); + int sv = e_src.v_start, ev = e_src.v_end; + // If this vertex belongs to a target face, use offset version + if (target_vertices.count(sv)) { + if (!normal_vremap.count(sv)) + normal_vremap[sv] = offset_vtx[sv]; + } else if (!normal_vremap.count(sv)) { + normal_vremap[sv] = result.add_vertex( + body.vertex_by_id(sv).point); + } + if (target_vertices.count(ev)) { + if (!normal_vremap.count(ev)) + normal_vremap[ev] = offset_vtx[ev]; + } else if (!normal_vremap.count(ev)) { + normal_vremap[ev] = result.add_vertex( + body.vertex_by_id(ev).point); + } + new_es.push_back(result.add_edge( + normal_vremap[sv], normal_vremap[ev])); + } + int new_lp = result.add_loop(new_es, lop.is_outer); + int new_surf = result.add_surface(body.surface(f_src.surface_id)); + result.add_face(new_surf, {new_lp}); + } + } + + // ── Fill corners: for vertices with 3+ target faces, add spherical patch ── + // Count target faces per vertex + std::map vertex_face_count; + for (int fi : face_ids) { + auto f_edges = body.face_edges(fi); + std::set seen; + for (int ei : f_edges) { + const auto& e = body.edge(ei); + if (!seen.count(e.v_start)) { + vertex_face_count[e.v_start]++; + seen.insert(e.v_start); + } + if (!seen.count(e.v_end)) { + vertex_face_count[e.v_end]++; + seen.insert(e.v_end); + } + } + } + + // Add spherical patches at corners touched by 3+ target faces + for (const auto& [vid, count] : vertex_face_count) { + if (count < 3) continue; + // Gather normals of target faces at this vertex + std::vector corner_normals; + for (int fi : face_ids) { + auto f_edges = body.face_edges(fi); + bool found = false; + for (int ei : f_edges) { + const auto& e = body.edge(ei); + if (e.v_start == vid || e.v_end == vid) { found = true; break; } + } + if (found) corner_normals.push_back(face_normal(body, fi)); + } + if (corner_normals.size() < 3) continue; + + Point3D orig = vertex_pos[vid]; + Point3D oc = { + orig.x() - vertex_normal[vid].x() * radius, + orig.y() - vertex_normal[vid].y() * radius, + orig.z() - vertex_normal[vid].z() * radius + }; + + // Build a triangular blend patch connecting the 3 offset edges + // using the first 3 normals for a spherical cap approximation + auto patch = make_sphere_corner_patch(oc, + corner_normals[0], corner_normals[1], corner_normals[2], radius); + + // Create 3 offset corner vertices + Point3D c0 = {oc.x() - corner_normals[0].x() * radius, + oc.y() - corner_normals[0].y() * radius, + oc.z() - corner_normals[0].z() * radius}; + Point3D c1 = {oc.x() - corner_normals[1].x() * radius, + oc.y() - corner_normals[1].y() * radius, + oc.z() - corner_normals[1].z() * radius}; + Point3D c2 = {oc.x() - corner_normals[2].x() * radius, + oc.y() - corner_normals[2].y() * radius, + oc.z() - corner_normals[2].z() * radius}; + + add_tri(result, c0, c1, c2, patch); + } + + finish_model(result, "multi_face_blend"); + return result; +} + +// ═══════════════════════════════════════════════ +// rolling_ball_blend +// ═══════════════════════════════════════════════ + +BrepModel rolling_ball_blend(const BrepModel& body, int edge_id, + double radius, int samples) { + int num_edges = static_cast(body.num_edges()); + if (edge_id < 0 || edge_id >= num_edges) return body; + if (radius <= 0.0) return body; + if (samples < 4) samples = 4; + + // 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; + + const auto& edge = body.edge(edge_id); + auto curve = edge_curve(body, edge_id); + + Vector3D n_a = face_normal(body, face_a); + Vector3D n_b = face_normal(body, face_b); + + // For rolling ball, compute the ball center trace + // Ball center = edge_point + r * bisector / cos(half_angle) + Vector3D bisector = (n_a + n_b).normalized(); + double cos_half = bisector.dot(n_a); + if (std::abs(cos_half) < 1e-8) return body; + double ball_offset = radius / cos_half; + + auto [t_min, t_max] = curve.domain(); + int N = samples; + + // Sample ball centers and contact points + std::vector centers(N + 1); + std::vector t1_pts(N + 1), t2_pts(N + 1); + + for (int i = 0; i <= N; ++i) { + double t_val = static_cast(i) / N; + double t = t_min + (t_max - t_min) * t_val; + Point3D ep = curve.evaluate(t); + + centers[i] = { + ep.x() + ball_offset * bisector.x(), + ep.y() + ball_offset * bisector.y(), + ep.z() + ball_offset * bisector.z() + }; + + // Contact points on each face (foot of perpendicular from center) + t1_pts[i] = { + centers[i].x() - radius * n_a.x(), + centers[i].y() - radius * n_a.y(), + centers[i].z() - radius * n_a.z() + }; + t2_pts[i] = { + centers[i].x() - radius * n_b.x(), + centers[i].y() - radius * n_b.y(), + centers[i].z() - radius * n_b.z() + }; + } + + BrepModel result; + + // ── Rebuild face_a: follow loop order ── + { + 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: follow loop order ── + { + 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 rolling ball envelope surface (piecewise patches) ── + for (int i = 0; i < N; ++i) { + // Envelope midpoint: on the sphere surface + Point3D mid_a = { + (t1_pts[i].x() + t1_pts[i+1].x()) * 0.5, + (t1_pts[i].y() + t1_pts[i+1].y()) * 0.5, + (t1_pts[i].z() + t1_pts[i+1].z()) * 0.5 + }; + Point3D mid_b = { + (t2_pts[i].x() + t2_pts[i+1].x()) * 0.5, + (t2_pts[i].y() + t2_pts[i+1].y()) * 0.5, + (t2_pts[i].z() + t2_pts[i+1].z()) * 0.5 + }; + // Arc peak: offset from center toward bisector + Point3D peak_a = { + centers[i].x() - radius * bisector.x(), + centers[i].y() - radius * bisector.y(), + centers[i].z() - radius * bisector.z() + }; + Point3D peak_b = { + centers[i+1].x() - radius * bisector.x(), + centers[i+1].y() - radius * bisector.y(), + centers[i+1].z() - radius * bisector.z() + }; + + auto patch = make_blend_patch(t1_pts[i], peak_a, t2_pts[i], + t1_pts[i+1], peak_b, t2_pts[i+1]); + add_quad(result, t1_pts[i], t1_pts[i+1], t2_pts[i+1], t2_pts[i], patch); + } + + finish_model(result, "rolling_ball_blend"); + return result; +} + +// ═══════════════════════════════════════════════ +// face_face_blend +// ═══════════════════════════════════════════════ + +BrepModel face_face_blend(const BrepModel& body, int face_a, int face_b, + double radius, int samples) { + int num_faces = static_cast(body.num_faces()); + if (face_a < 0 || face_a >= num_faces) return body; + if (face_b < 0 || face_b >= num_faces) return body; + if (radius <= 0.0) return body; + if (samples < 2) samples = 2; + + // Find shared edges between face_a and face_b + auto edges_a = body.face_edges(face_a); + auto edges_b = body.face_edges(face_b); + + std::vector shared_edges; + for (int ea : edges_a) { + for (int eb : edges_b) { + if (edges_same(body, ea, eb)) { + shared_edges.push_back(ea); + break; + } + } + } + + if (shared_edges.empty()) { + // No shared edge — fall back to single-edge blend on first common edge + // Try edge_faces to find an edge touching both + for (int ea : edges_a) { + auto efs = body.edge_faces(ea); + bool touches_b = false; + for (int f : efs) { + if (f == face_b) { touches_b = true; break; } + } + if (touches_b) { shared_edges.push_back(ea); break; } + } + } + + if (shared_edges.empty()) return body; + + // For simplicity, blend along the first shared edge + int edge_id = shared_edges[0]; + const auto& edge = body.edge(edge_id); + auto curve = edge_curve(body, edge_id); + + 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; + double offset_dist = radius / cos_half; + + auto [t_min, t_max] = curve.domain(); + int N = samples; + + // Sample along the shared edge + 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 t = t_min + (t_max - t_min) * t_val; + Point3D ep = curve.evaluate(t); + + Point3D C = { + ep.x() + offset_dist * bisector.x(), + ep.y() + offset_dist * bisector.y(), + ep.z() + offset_dist * bisector.z() + }; + + // Tangent points on each face + t1_pts[i] = { + C.x() - radius * n_a.x(), C.y() - radius * n_a.y(), C.z() - radius * n_a.z() + }; + t2_pts[i] = { + C.x() - radius * n_b.x(), C.y() - radius * n_b.y(), C.z() - radius * n_b.z() + }; + mid_pts[i] = C; + } + + BrepModel result; + + // ── Rebuild face_a: follow loop order ── + { + 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: follow loop order ── + { + 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 constant-radius blend surface strips ── + 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); + } + + finish_model(result, "face_face_blend"); + return result; +} + +} // namespace vde::brep diff --git a/src/curves/surface_analysis.cpp b/src/curves/surface_analysis.cpp new file mode 100644 index 0000000..161de92 --- /dev/null +++ b/src/curves/surface_analysis.cpp @@ -0,0 +1,537 @@ +#include "vde/curves/surface_analysis.h" +#include +#include +#include + +namespace vde::curves { + +using core::Point3D; +using core::Vector3D; + +// ═══════════════════════════════════════════════════════════ +// Internal Helpers +// ═══════════════════════════════════════════════════════════ + +namespace { + +/// Compute the first and second fundamental form coefficients at a surface point +/// I = [E F; F G], II = [L M; M N] +struct SurfaceDerivs { + Vector3D Su, Sv; // first derivatives + Vector3D Suu, Suv, Svv; // second derivatives + Vector3D N; // unit normal + double E, F, G; // first fundamental form + double L, M, Ncoeff; // second fundamental form +}; + +SurfaceDerivs compute_surface_derivs(const NurbsSurface& surf, double u, double v, + double eps = 1e-5) { + SurfaceDerivs d; + + d.Su = surf.derivative_u(u, v); + d.Sv = surf.derivative_v(u, v); + d.N = surf.normal(u, v); + + // Second derivatives via finite difference + d.Suu = (surf.derivative_u(u + eps, v) - surf.derivative_u(u - eps, v)) / (2.0 * eps); + d.Suv = (surf.derivative_u(u, v + eps) - surf.derivative_u(u, v - eps)) / (2.0 * eps); + d.Svv = (surf.derivative_v(u, v + eps) - surf.derivative_v(u, v - eps)) / (2.0 * eps); + + // First fundamental form + d.E = d.Su.dot(d.Su); + d.F = d.Su.dot(d.Sv); + d.G = d.Sv.dot(d.Sv); + + // Second fundamental form + d.L = d.Suu.dot(d.N); + d.M = d.Suv.dot(d.N); + d.Ncoeff = d.Svv.dot(d.N); + + return d; +} + +/// Compute principal curvatures from fundamental forms +/// κ satisfies: (EG-F²)κ² - (EN+GL-2FM)κ + (LN-M²) = 0 +void principal_curvatures(const SurfaceDerivs& d, double& k1, double& k2, + Vector3D& dir1, Vector3D& dir2) { + double detI = d.E * d.G - d.F * d.F; + if (std::abs(detI) < 1e-16) { + k1 = k2 = 0.0; + dir1 = d.Su.normalized(); + dir2 = d.Sv.normalized(); + return; + } + + // Mean curvature H, Gaussian curvature K + double H = (d.E * d.Ncoeff - 2.0 * d.F * d.M + d.G * d.L) / (2.0 * detI); + double K = (d.L * d.Ncoeff - d.M * d.M) / detI; + + double disc = H * H - K; + if (disc < 0) disc = 0; + double sqrt_disc = std::sqrt(disc); + + k1 = H + sqrt_disc; + k2 = H - sqrt_disc; + + // Principal directions (in parameter space) + // Direction for k1: solve (L - k1*E) du + (M - k1*F) dv = 0 + auto principal_dir = [&](double k) -> Vector3D { + double a = d.L - k * d.E; + double b = d.M - k * d.F; + // If a≈0 and b≈0, any direction is principal (umbilic point) + if (std::abs(a) < 1e-12 && std::abs(b) < 1e-12) { + return d.Su.normalized(); + } + // Find direction in tangent plane + if (std::abs(a) > std::abs(b)) { + double dv_du = -a / b; + Vector3D dir = d.Su + dv_du * d.Sv; + return dir.normalized(); + } else { + double du_dv = -b / a; + Vector3D dir = du_dv * d.Su + d.Sv; + return dir.normalized(); + } + }; + + dir1 = principal_dir(k1); + dir2 = principal_dir(k2); + + // Ensure dir1, dir2 are orthogonal in tangent plane + if (std::abs(dir1.dot(dir2)) > 0.9) { + // Near-umbilic: use Su, Sv as approximate principal directions + dir1 = d.Su.normalized(); + // Gram-Schmidt for dir2 + dir2 = (d.Sv - d.Sv.dot(dir1) * dir1).normalized(); + } +} + +/// Euclidean distance between two 3D points +inline double dist3d(const Point3D& a, const Point3D& b) { + return (a - b).norm(); +} + +} // namespace + +// ═══════════════════════════════════════════════════════════ +// CurvatureMap +// ═══════════════════════════════════════════════════════════ + +bool CurvatureMap::is_developable(double tol) const { + return max_gaussian < tol && min_gaussian > -tol; +} + +CurvatureMap curvature_map(const NurbsSurface& surf, int res_u, int res_v) { + CurvatureMap 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]; + + result.grid.resize(res_u + 1); + for (int i = 0; i <= res_u; ++i) { + result.grid[i].resize(res_v + 1); + 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 derivs = compute_surface_derivs(surf, u, v); + CurvaturePoint& cp = result.grid[i][j]; + + principal_curvatures(derivs, cp.k1, cp.k2, cp.d1, cp.d2); + cp.gaussian = cp.k1 * cp.k2; + cp.mean = 0.5 * (cp.k1 + cp.k2); + + result.min_gaussian = std::min(result.min_gaussian, cp.gaussian); + result.max_gaussian = std::max(result.max_gaussian, cp.gaussian); + result.min_mean = std::min(result.min_mean, cp.mean); + result.max_mean = std::max(result.max_mean, cp.mean); + } + } + + return result; +} + +// ═══════════════════════════════════════════════════════════ +// Zebra Stripe +// ═══════════════════════════════════════════════════════════ + +ZebraStripeResult zebra_stripe( + const NurbsSurface& surf, + const Vector3D& light_dir, + int res_u, int res_v) { + + ZebraStripeResult result; + result.res_u = res_u; + result.res_v = res_v; + + Vector3D L = light_dir.normalized(); + result.light_direction = L; + + 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)); + + 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); + double dot = N.dot(L); + + // Zebra stripe: sinusoidal pattern based on elevation angle + // Maps [-1,1] dot range to [0,1] intensity using smoothstep-like pattern + double intensity = 0.5 * (1.0 + std::sin(dot * 20.0 * M_PI)); + result.intensity[i][j] = std::max(0.0, std::min(1.0, intensity)); + } + } + + return result; +} + +// ═══════════════════════════════════════════════════════════ +// Curvature Comb +// ═══════════════════════════════════════════════════════════ + +CurvatureComb curvature_comb(const NurbsCurve& curve, int res, double scale) { + CurvatureComb result; + + auto dom = curve.domain(); + double t_min = dom.first; + double t_max = dom.second; + + if (t_max <= t_min || res < 2) return result; + + // First pass: find curvature range for auto-scaling + for (int i = 0; i <= res; ++i) { + double t = t_min + (t_max - t_min) * i / res; + Vector3D d1 = curve.derivative(t, 1); + Vector3D d2 = curve.derivative(t, 2); + + double speed2 = d1.squaredNorm(); + if (speed2 < 1e-18) continue; + + // Curvature magnitude: κ = |r' × r''| / |r'|³ + double k = d1.cross(d2).norm() / (speed2 * std::sqrt(speed2)); + + result.min_curvature = std::min(result.min_curvature, k); + result.max_curvature = std::max(result.max_curvature, k); + } + + // Auto-scale: make max comb length = some fraction of curve extent + if (scale <= 0.0 && result.max_curvature > 1e-12) { + // Estimate curve extent + double extent = (curve.evaluate(t_max) - curve.evaluate(t_min)).norm(); + scale = extent * 0.15 / result.max_curvature; + } + if (scale <= 0.0) scale = 1.0; + result.scale = scale; + + // Second pass: compute comb points + for (int i = 0; i <= res; ++i) { + double t = t_min + (t_max - t_min) * i / res; + + Point3D pt = curve.evaluate(t); + Vector3D d1 = curve.derivative(t, 1); + Vector3D d2 = curve.derivative(t, 2); + + double speed2 = d1.squaredNorm(); + Vector3D curve_normal; + + if (speed2 < 1e-18) { + curve_normal = Vector3D::Zero(); + } else { + // Frenet normal: N = (r' × r'') × r' / (|r'| · |r' × r''|) + Vector3D cross = d1.cross(d2); + double cross_norm = cross.norm(); + if (cross_norm < 1e-18) { + curve_normal = Vector3D::Zero(); + } else { + curve_normal = cross.cross(d1).normalized(); + } + } + + double k = 0.0; + if (speed2 > 1e-18) { + k = d1.cross(d2).norm() / (speed2 * std::sqrt(speed2)); + } + + CurvatureCombPoint cp; + cp.curve_point = pt; + cp.curvature = k; + cp.curvature_vec = curve_normal * k * scale; + + result.points.push_back(cp); + } + + return result; +} + +// ═══════════════════════════════════════════════════════════ +// Deviation Analysis +// ═══════════════════════════════════════════════════════════ + +DeviationResult deviation_analysis( + const NurbsSurface& surf_a, + const NurbsSurface& surf_b, + int samples) { + + DeviationResult result; + + const auto& ku_a = surf_a.knots_u(); + const auto& kv_a = surf_a.knots_v(); + + if (ku_a.empty() || kv_a.empty() || samples < 2) return result; + + int pu = surf_a.degree_u(); + int pv = surf_a.degree_v(); + + double u_min = ku_a[pu]; + double u_max = ku_a[ku_a.size() - 1 - pu]; + double v_min = kv_a[pv]; + double v_max = kv_a[kv_a.size() - 1 - pv]; + + result.sample_count = static_cast(samples * samples); + + double sum_sq = 0.0; + double sum_abs = 0.0; + + for (int i = 0; i < samples; ++i) { + double u = u_min + (u_max - u_min) * i / (samples - 1); + for (int j = 0; j < samples; ++j) { + double v = v_min + (v_max - v_min) * j / (samples - 1); + + Point3D pa = surf_a.evaluate(u, v); + Vector3D na = surf_a.normal(u, v); + + // Find closest point on surf_b by sampling surf_b parameter domain + const auto& ku_b = surf_b.knots_u(); + const auto& kv_b = surf_b.knots_v(); + if (ku_b.empty() || kv_b.empty()) continue; + + int pu_b = surf_b.degree_u(); + int pv_b = surf_b.degree_v(); + double ub_min = ku_b[pu_b]; + double ub_max = ku_b[ku_b.size() - 1 - pu_b]; + double vb_min = kv_b[pv_b]; + double vb_max = kv_b[kv_b.size() - 1 - pv_b]; + + // Coarse-to-fine closest point search on surf_b + int search_res = 20; + double best_dist = std::numeric_limits::max(); + double best_ub = ub_min, best_vb = vb_min; + + for (int si = 0; si <= search_res; ++si) { + double ub = ub_min + (ub_max - ub_min) * si / search_res; + for (int sj = 0; sj <= search_res; ++sj) { + double vb = vb_min + (vb_max - vb_min) * sj / search_res; + Point3D pb = surf_b.evaluate(ub, vb); + double d = (pa - pb).norm(); + if (d < best_dist) { + best_dist = d; + best_ub = ub; + best_vb = vb; + } + } + } + + Point3D pb = surf_b.evaluate(best_ub, best_vb); + + // Signed deviation: positive if pb is on the positive normal side of pa + Vector3D diff = pb - pa; + double dev = diff.dot(na); + + result.deviations.push_back(dev); + result.sample_points_a.push_back(pa); + result.sample_points_b.push_back(pb); + + if (dev > 0) { + result.max_positive = std::max(result.max_positive, dev); + } else { + result.max_negative = std::min(result.max_negative, dev); + } + + sum_sq += dev * dev; + sum_abs += std::abs(dev); + } + } + + size_t n = result.deviations.size(); + if (n > 0) { + result.rms = std::sqrt(sum_sq / n); + result.mean_absolute = sum_abs / n; + } + + return result; +} + +// ═══════════════════════════════════════════════════════════ +// Draft Face Angle (Enhanced) +// ═══════════════════════════════════════════════════════════ + +double draft_face_angle( + const NurbsSurface& surf, + double u, double v, + const Vector3D& pull_dir) { + + Vector3D dir = pull_dir.normalized(); + Vector3D N = surf.normal(u, v); + double dot = N.dot(dir); + + // Draft angle = π/2 - arccos(|n·d|) + // Positive: face opens outward along pull direction + double base_angle = M_PI_2 - std::acos(std::max(-1.0, std::min(1.0, std::abs(dot)))); + return (dot < 0) ? base_angle : -base_angle; +} + +double draft_face_angle( + const brep::BrepModel& body, + int face_id, + const Vector3D& pull_dir, + int samples) { + + if (face_id < 0 || face_id >= static_cast(body.num_faces())) { + return 0.0; + } + + const auto& face = body.face(face_id); + if (face.surface_id < 0 || face.surface_id >= static_cast(body.num_surfaces())) { + return 0.0; + } + + const auto& surf = body.surface(face.surface_id); + const auto& ku = surf.knots_u(); + const auto& kv = surf.knots_v(); + + if (ku.empty() || kv.empty()) return 0.0; + + 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]; + + // Area-weighted sampling + double total_weight = 0.0; + double weighted_sum = 0.0; + + for (int i = 0; i < samples; ++i) { + double u = u_min + (u_max - u_min) * i / (samples - 1); + for (int j = 0; j < samples; ++j) { + double v = v_min + (v_max - v_min) * j / (samples - 1); + + // Area element weight: |Su × Sv| + Vector3D Su = surf.derivative_u(u, v); + Vector3D Sv = surf.derivative_v(u, v); + double area_weight = Su.cross(Sv).norm(); + + double angle = draft_face_angle(surf, u, v, pull_dir); + + // Apply face reversal + if (face.reversed) angle = -angle; + + weighted_sum += angle * area_weight; + total_weight += area_weight; + } + } + + if (total_weight > 1e-15) { + return weighted_sum / total_weight; + } + return 0.0; +} + +DraftAngleDistribution draft_face_angle_distribution( + const brep::BrepModel& body, + int face_id, + const Vector3D& pull_dir, + int samples) { + + DraftAngleDistribution result; + result.min_angle = std::numeric_limits::max(); + result.max_angle = -std::numeric_limits::max(); + + if (face_id < 0 || face_id >= static_cast(body.num_faces())) return result; + + const auto& face = body.face(face_id); + if (face.surface_id < 0 || face.surface_id >= static_cast(body.num_surfaces())) return result; + + const auto& surf = body.surface(face.surface_id); + 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]; + + double sum_angle = 0.0; + double total_weight = 0.0; + double weighted_sum = 0.0; + + for (int i = 0; i < samples; ++i) { + double u = u_min + (u_max - u_min) * i / (samples - 1); + for (int j = 0; j < samples; ++j) { + double v = v_min + (v_max - v_min) * j / (samples - 1); + + double angle = draft_face_angle(surf, u, v, pull_dir); + if (face.reversed) angle = -angle; + + Vector3D Su = surf.derivative_u(u, v); + Vector3D Sv = surf.derivative_v(u, v); + double area_weight = Su.cross(Sv).norm(); + + result.sample_angles.push_back(angle); + result.sample_params.push_back({u, v}); + + sum_angle += angle; + weighted_sum += angle * area_weight; + total_weight += area_weight; + + result.min_angle = std::min(result.min_angle, angle); + result.max_angle = std::max(result.max_angle, angle); + } + } + + size_t n = result.sample_angles.size(); + result.mean_angle = (n > 0) ? sum_angle / n : 0.0; + result.area_weighted_angle = (total_weight > 1e-15) ? weighted_sum / total_weight : 0.0; + + if (result.sample_angles.empty()) { + result.min_angle = 0.0; + result.max_angle = 0.0; + } + + return result; +} + +} // namespace vde::curves diff --git a/src/curves/surface_continuity.cpp b/src/curves/surface_continuity.cpp new file mode 100644 index 0000000..b7d016b --- /dev/null +++ b/src/curves/surface_continuity.cpp @@ -0,0 +1,610 @@ +#include "vde/curves/surface_continuity.h" +#include +#include +#include + +namespace vde::curves { + +using core::Point3D; +using core::Vector3D; + +// ═══════════════════════════════════════════════════════════ +// Utility +// ═══════════════════════════════════════════════════════════ + +const char* to_string(ContinuityLevel level) { + switch (level) { + case ContinuityLevel::None: return "None"; + case ContinuityLevel::G0: return "G0"; + case ContinuityLevel::G1: return "G1"; + case ContinuityLevel::G2: return "G2"; + case ContinuityLevel::G3: return "G3"; + } + return "Unknown"; +} + +namespace { + +/// Check if vector is zero (within tolerance) +bool is_zero_vector(const Vector3D& v, double tol = 1e-12) { + return v.norm() < tol; +} + +/// Check if two vectors are collinear and point in the same direction +bool are_collinear_same_direction(const Vector3D& a, const Vector3D& b, double tol = 1e-9) { + if (is_zero_vector(a) || is_zero_vector(b)) return false; + double cross_norm = a.cross(b).norm(); + double dot = a.dot(b); + return cross_norm < tol && dot > 0; +} + +/// Compute curvature vector: κ = (r' × r'') × r' / |r'|^4 +/// Returns the curvature vector (magnitude = curvature, direction = normal) +Vector3D curvature_vector(const Vector3D& d1, const Vector3D& d2) { + double speed2 = d1.squaredNorm(); + if (speed2 < 1e-18) return Vector3D::Zero(); + double speed4 = speed2 * speed2; + return (d1.cross(d2)).cross(d1) / speed4; +} + +/// Compute torsion (rate of change of binormal along the curve) +/// τ = (r' × r'') · r''' / |r' × r''|² +double torsion(const Vector3D& d1, const Vector3D& d2, const Vector3D& d3) { + Vector3D cross12 = d1.cross(d2); + double cross_norm2 = cross12.squaredNorm(); + if (cross_norm2 < 1e-18) return 0.0; + return cross12.dot(d3) / cross_norm2; +} + +/// Check curvature continuity (G2): curvature vectors match +bool curvature_match(const Vector3D& d1a, const Vector3D& d2a, + const Vector3D& d1b, const Vector3D& d2b, double tol = 1e-6) { + auto kappa_a = curvature_vector(d1a, d2a); + auto kappa_b = curvature_vector(d1b, d2b); + return (kappa_a - kappa_b).norm() < tol; +} + +/// Check curvature variation continuity (G3): +/// Curvature derivative vectors match. +/// dκ/dt = (r'×(r'×r'''))×r' / |r'|^5 (simplified) +Vector3D curvature_derivative_vector(const Vector3D& d1, const Vector3D& d2, + const Vector3D& d3) { + double speed2 = d1.squaredNorm(); + if (speed2 < 1e-18) return Vector3D::Zero(); + double speed5 = std::pow(speed2, 2.5); + auto cross13 = d1.cross(d3); + auto outer = cross13.cross(d1); + return outer / speed5; +} + +bool curvature_derivative_match(const Vector3D& d1a, const Vector3D& d2a, + const Vector3D& d3a, + const Vector3D& d1b, const Vector3D& d2b, + const Vector3D& d3b, double tol = 1e-6) { + auto dkappa_a = curvature_derivative_vector(d1a, d2a, d3a); + auto dkappa_b = curvature_derivative_vector(d1b, d2b, d3b); + return (dkappa_a - dkappa_b).norm() < tol; +} + +} // namespace + +// ═══════════════════════════════════════════════════════════ +// Curve Continuity +// ═══════════════════════════════════════════════════════════ + +ContinuityLevel continuity_type( + const NurbsCurve& curve_a, double t_a, + const NurbsCurve& curve_b, double t_b, + double tol) { + + // --- G0: position --- + Point3D pa = curve_a.evaluate(t_a); + Point3D pb = curve_b.evaluate(t_b); + double dist = (pa - pb).norm(); + if (dist > tol) return ContinuityLevel::None; + + // --- G1: tangent collinear --- + Vector3D d1a = curve_a.derivative(t_a, 1); + Vector3D d1b = curve_b.derivative(t_b, 1); + + // If either derivative is degenerate, can't determine G1+ + if (is_zero_vector(d1a) || is_zero_vector(d1b)) return ContinuityLevel::G0; + + if (!are_collinear_same_direction(d1a, d1b, tol)) { + return ContinuityLevel::G0; + } + + // --- G2: curvature continuity --- + Vector3D d2a = curve_a.derivative(t_a, 2); + Vector3D d2b = curve_b.derivative(t_b, 2); + + if (!curvature_match(d1a, d2a, d1b, d2b, tol)) { + return ContinuityLevel::G1; + } + + // --- G3: curvature variation continuity --- + Vector3D d3a = curve_a.derivative(t_a, 3); + Vector3D d3b = curve_b.derivative(t_b, 3); + + if (!curvature_derivative_match(d1a, d2a, d3a, d1b, d2b, d3b, tol)) { + return ContinuityLevel::G2; + } + + return ContinuityLevel::G3; +} + +ContinuityLevel continuity_type( + const NurbsCurve& curve_a, + const NurbsCurve& curve_b, + double tol) { + + auto dom_a = curve_a.domain(); + auto dom_b = curve_b.domain(); + + Point3D a_start = curve_a.evaluate(dom_a.first); + Point3D a_end = curve_a.evaluate(dom_a.second); + Point3D b_start = curve_b.evaluate(dom_b.first); + Point3D b_end = curve_b.evaluate(dom_b.second); + + // Try all endpoint pairs + struct Candidate { double ta, tb; }; + std::vector candidates; + + auto local_tol = tol * 10.0; // relaxed for endpoint matching + if ((a_end - b_start).norm() < local_tol) + candidates.push_back({dom_a.second, dom_b.first}); + if ((a_end - b_end).norm() < local_tol) + candidates.push_back({dom_a.second, dom_b.second}); + if ((a_start - b_start).norm() < local_tol) + candidates.push_back({dom_a.first, dom_b.first}); + if ((a_start - b_end).norm() < local_tol) + candidates.push_back({dom_a.first, dom_b.second}); + + // If no endpoints match, try the other way (b→a) + if (candidates.empty()) { + return ContinuityLevel::None; + } + + // Return the best continuity found among all matching endpoint pairs + ContinuityLevel best = ContinuityLevel::None; + for (const auto& c : candidates) { + auto level = continuity_type(curve_a, c.ta, curve_b, c.tb, tol); + if (static_cast(level) > static_cast(best)) { + best = level; + } + } + return best; +} + +// ═══════════════════════════════════════════════════════════ +// Surface Continuity +// ═══════════════════════════════════════════════════════════ + +ContinuityLevel surface_continuity( + const NurbsSurface& surf_a, + const NurbsSurface& surf_b, + std::pair, std::pair> edge_param_range_a, + std::pair, std::pair> edge_param_range_b, + int samples, + double tol) { + + double ua_start = edge_param_range_a.first.first; + double va_start = edge_param_range_a.first.second; + double ua_end = edge_param_range_a.second.first; + double va_end = edge_param_range_a.second.second; + + double ub_start = edge_param_range_b.first.first; + double vb_start = edge_param_range_b.first.second; + double ub_end = edge_param_range_b.second.first; + double vb_end = edge_param_range_b.second.second; + + ContinuityLevel worst = ContinuityLevel::G3; + + for (int i = 0; i <= samples; ++i) { + double t = static_cast(i) / samples; + + double ua = ua_start + t * (ua_end - ua_start); + double va = va_start + t * (va_end - va_start); + double ub = ub_start + t * (ub_end - ub_start); + double vb = vb_start + t * (vb_end - vb_start); + + // G0: position + Point3D pa = surf_a.evaluate(ua, va); + Point3D pb = surf_b.evaluate(ub, vb); + if ((pa - pb).norm() > tol) { + if (static_cast(ContinuityLevel::None) < static_cast(worst)) + worst = ContinuityLevel::None; + continue; + } + + // Determine which direction is "along the edge" to compare derivatives + // If u varies along edge (i.e. ua_end != ua_start), derivatives are du + // If v varies along edge, derivatives are dv + bool along_u_a = std::abs(ua_end - ua_start) > std::abs(va_end - va_start); + bool along_u_b = std::abs(ub_end - ub_start) > std::abs(vb_end - vb_start); + + Vector3D d1a, d1b; + if (along_u_a) d1a = surf_a.derivative_u(ua, va); + else d1a = surf_a.derivative_v(ua, va); + if (along_u_b) d1b = surf_b.derivative_u(ub, vb); + else d1b = surf_b.derivative_v(ub, vb); + + if (is_zero_vector(d1a) || is_zero_vector(d1b)) { + if (static_cast(ContinuityLevel::G0) < static_cast(worst)) + worst = ContinuityLevel::G0; + continue; + } + + if (!are_collinear_same_direction(d1a, d1b, tol)) { + if (static_cast(ContinuityLevel::G0) < static_cast(worst)) + worst = ContinuityLevel::G0; + continue; + } + + // G2: second derivative + // For surfaces, we need the second derivative along the edge direction + // We approximate using finite differences along the edge + double eps = 1e-4; + double ua_mid = ua, va_mid = va; + double ub_mid = ub, vb_mid = vb; + + // Compute second derivatives via finite difference along the edge direction + auto second_deriv_along = [&](const NurbsSurface& s, double u, double v, + bool along_u, double du_step, double dv_step) -> Vector3D { + auto fwd = s.evaluate(u + du_step, v + dv_step); + auto mid = s.evaluate(u, v); + auto bwd = s.evaluate(u - du_step, v - dv_step); + double step_sq = (du_step * du_step + dv_step * dv_step); + if (step_sq < 1e-20) return Vector3D::Zero(); + return (fwd - 2.0 * mid + bwd) / step_sq; + }; + + double du_step_a = along_u_a ? eps : 0.0; + double dv_step_a = along_u_a ? 0.0 : eps; + double du_step_b = along_u_b ? eps : 0.0; + double dv_step_b = along_u_b ? 0.0 : eps; + + // Use actual derivative for d2 + Vector3D d2a, d2b; + if (along_u_a) + d2a = surf_a.derivative_u(ua, va); // This is du, we need duu — use finite diff + else + d2a = surf_a.derivative_v(ua, va); // This is dv, we need dvv — use finite diff + + // Actually, the proper approach: get the actual second derivative + // We'll use the derivative methods and finite difference approximation + // For NurbsSurface derivative_u/v gives first derivatives + // We need second derivatives along the edge direction + // Use finite differences for second derivative + { + double h = eps; + Vector3D d1_fwd, d1_bwd; + if (along_u_a) { + d1_fwd = surf_a.derivative_u(ua + h, va); + d1_bwd = surf_a.derivative_u(ua - h, va); + } else { + d1_fwd = surf_a.derivative_v(ua, va + h); + d1_bwd = surf_a.derivative_v(ua, va - h); + } + d2a = (d1_fwd - d1_bwd) / (2.0 * h); + } + { + double h = eps; + Vector3D d1_fwd, d1_bwd; + if (along_u_b) { + d1_fwd = surf_b.derivative_u(ub + h, vb); + d1_bwd = surf_b.derivative_u(ub - h, vb); + } else { + d1_fwd = surf_b.derivative_v(ub, vb + h); + d1_bwd = surf_b.derivative_v(ub, vb - h); + } + d2b = (d1_fwd - d1_bwd) / (2.0 * h); + } + + if (!curvature_match(d1a, d2a, d1b, d2b, tol)) { + if (static_cast(ContinuityLevel::G1) < static_cast(worst)) + worst = ContinuityLevel::G1; + continue; + } + + // G3: third derivative (finite difference of second derivative) + Vector3D d3a, d3b; + { + double h = eps * 2.0; + Vector3D d2_fwd, d2_bwd; + if (along_u_a) { + // compute second deriv at ua+h and ua-h + auto d1_p2 = surf_a.derivative_u(ua + 2*h, va); + auto d1_0 = surf_a.derivative_u(ua, va); + auto d1_n2 = surf_a.derivative_u(ua - 2*h, va); + d2_fwd = (d1_p2 - d1_0) / (2.0 * h); + d2_bwd = (d1_0 - d1_n2) / (2.0 * h); + } else { + auto d1_p2 = surf_a.derivative_v(ua, va + 2*h); + auto d1_0 = surf_a.derivative_v(ua, va); + auto d1_n2 = surf_a.derivative_v(ua, va - 2*h); + d2_fwd = (d1_p2 - d1_0) / (2.0 * h); + d2_bwd = (d1_0 - d1_n2) / (2.0 * h); + } + d3a = (d2_fwd - d2_bwd) / (2.0 * h); + } + { + double h = eps * 2.0; + Vector3D d2_fwd, d2_bwd; + if (along_u_b) { + auto d1_p2 = surf_b.derivative_u(ub + 2*h, vb); + auto d1_0 = surf_b.derivative_u(ub, vb); + auto d1_n2 = surf_b.derivative_u(ub - 2*h, vb); + d2_fwd = (d1_p2 - d1_0) / (2.0 * h); + d2_bwd = (d1_0 - d1_n2) / (2.0 * h); + } else { + auto d1_p2 = surf_b.derivative_v(ub, vb + 2*h); + auto d1_0 = surf_b.derivative_v(ub, vb); + auto d1_n2 = surf_b.derivative_v(ub, vb - 2*h); + d2_fwd = (d1_p2 - d1_0) / (2.0 * h); + d2_bwd = (d1_0 - d1_n2) / (2.0 * h); + } + d3b = (d2_fwd - d2_bwd) / (2.0 * h); + } + + if (!curvature_derivative_match(d1a, d2a, d3a, d1b, d2b, d3b, tol)) { + if (static_cast(ContinuityLevel::G2) < static_cast(worst)) + worst = ContinuityLevel::G2; + continue; + } + } + + return worst; +} + +ContinuityLevel surface_continuity( + const NurbsSurface& surf_a, + const NurbsSurface& surf_b, + int samples, + double tol) { + + const auto& ku_a = surf_a.knots_u(); + const auto& kv_a = surf_a.knots_v(); + const auto& ku_b = surf_b.knots_u(); + const auto& kv_b = surf_b.knots_v(); + + if (ku_a.empty() || kv_a.empty() || ku_b.empty() || kv_b.empty()) + return ContinuityLevel::None; + + int pu_a = surf_a.degree_u(), pv_a = surf_a.degree_v(); + int pu_b = surf_b.degree_u(), pv_b = surf_b.degree_v(); + + double ua_min = ku_a[pu_a], ua_max = ku_a[ku_a.size() - 1 - pu_a]; + double va_min = kv_a[pv_a], va_max = kv_a[kv_a.size() - 1 - pv_a]; + double ub_min = ku_b[pu_b], ub_max = ku_b[ku_b.size() - 1 - pu_b]; + double vb_min = kv_b[pv_b], vb_max = kv_b[kv_b.size() - 1 - pv_b]; + + struct Edge { + double u0, v0, u1, v1; + bool along_u; + }; + + std::vector edges_a = { + {ua_min, va_min, ua_min, va_max, false}, // u=umin edge + {ua_min, va_max, ua_max, va_max, true}, // v=vmax edge + {ua_max, va_max, ua_max, va_min, false}, // u=umax edge + {ua_max, va_min, ua_min, va_min, true}, // v=vmin edge + }; + + std::vector edges_b = { + {ub_min, vb_min, ub_min, vb_max, false}, + {ub_min, vb_max, ub_max, vb_max, true}, + {ub_max, vb_max, ub_max, vb_min, false}, + {ub_max, vb_min, ub_min, vb_min, true}, + }; + + ContinuityLevel best = ContinuityLevel::None; + + // Compare edges and align parameterization direction + double relaxed_tol = tol * 100.0; + for (const auto& ea : edges_a) { + Point3D ma_start = surf_a.evaluate(ea.u0, ea.v0); + Point3D ma_end = surf_a.evaluate(ea.u1, ea.v1); + for (const auto& eb : edges_b) { + Point3D mb_start = surf_b.evaluate(eb.u0, eb.v0); + Point3D mb_end = surf_b.evaluate(eb.u1, eb.v1); + + // Check if edges match (midpoint comparison) + Point3D ma_mid = surf_a.evaluate(0.5*(ea.u0+ea.u1), 0.5*(ea.v0+ea.v1)); + Point3D mb_mid = surf_b.evaluate(0.5*(eb.u0+eb.u1), 0.5*(eb.v0+eb.v1)); + if ((ma_mid - mb_mid).norm() >= relaxed_tol) + continue; + + // Align edge directions: if start→start is farther than start→end, reverse b + double d_ss = (ma_start - mb_start).norm(); + double d_se = (ma_start - mb_end).norm(); + std::pair,std::pair> eb_range; + + if (d_se < d_ss) { + // Reverse b: swap start/end + eb_range = {{eb.u1, eb.v1}, {eb.u0, eb.v0}}; + } else { + eb_range = {{eb.u0, eb.v0}, {eb.u1, eb.v1}}; + } + + auto level = surface_continuity( + surf_a, surf_b, + {{ea.u0, ea.v0}, {ea.u1, ea.v1}}, + eb_range, + samples, tol); + if (static_cast(level) > static_cast(best)) + best = level; + } + } + + return best; +} + +SurfaceContinuityReport surface_continuity_report( + const NurbsSurface& surf_a, + const NurbsSurface& surf_b, + int samples, + double tol) { + + SurfaceContinuityReport report; + report.worst = ContinuityLevel::G3; + report.best = ContinuityLevel::None; + report.g1_ratio = 0.0; + report.g2_ratio = 0.0; + + const auto& ku_a = surf_a.knots_u(); + const auto& kv_a = surf_a.knots_v(); + const auto& ku_b = surf_b.knots_u(); + const auto& kv_b = surf_b.knots_v(); + + if (ku_a.empty() || kv_a.empty() || ku_b.empty() || kv_b.empty()) + return report; + + int pu_a = surf_a.degree_u(), pv_a = surf_a.degree_v(); + int pu_b = surf_b.degree_u(), pv_b = surf_b.degree_v(); + + double ua_min = ku_a[pu_a], ua_max = ku_a[ku_a.size() - 1 - pu_a]; + double va_min = kv_a[pv_a], va_max = kv_a[kv_a.size() - 1 - pv_a]; + double ub_min = ku_b[pu_b], ub_max = ku_b[ku_b.size() - 1 - pu_b]; + double vb_min = kv_b[pv_b], vb_max = kv_b[kv_b.size() - 1 - pv_b]; + + // Find best-matching edge pair + struct EdgeInfo { + double u0, v0, u1, v1; + }; + + // Compare all edge pairs, pick the best matching one + std::vector edges_a = { + {ua_min, va_min, ua_min, va_max}, + {ua_min, va_max, ua_max, va_max}, + {ua_max, va_max, ua_max, va_min}, + {ua_max, va_min, ua_min, va_min}, + }; + std::vector edges_b = { + {ub_min, vb_min, ub_min, vb_max}, + {ub_min, vb_max, ub_max, vb_max}, + {ub_max, vb_max, ub_max, vb_min}, + {ub_max, vb_min, ub_min, vb_min}, + }; + + double best_dist = std::numeric_limits::max(); + std::pair best_pair; + bool found = false; + + double relaxed_tol = tol * 100.0; + for (const auto& ea : edges_a) { + Point3D ma = surf_a.evaluate(0.5*(ea.u0+ea.u1), 0.5*(ea.v0+ea.v1)); + for (const auto& eb : edges_b) { + Point3D mb = surf_b.evaluate(0.5*(eb.u0+eb.u1), 0.5*(eb.v0+eb.v1)); + double d = (ma - mb).norm(); + if (d < relaxed_tol && d < best_dist) { + best_dist = d; + best_pair = {ea, eb}; + found = true; + } + } + } + + if (!found) return report; + + const auto& ea = best_pair.first; + const auto& eb = best_pair.second; + + bool along_u_a = std::abs(ea.u1 - ea.u0) > std::abs(ea.v1 - ea.v0); + bool along_u_b = std::abs(eb.u1 - eb.u0) > std::abs(eb.v1 - eb.v0); + double eps = 1e-4; + + int g1_count = 0, g2_count = 0; + + for (int i = 0; i <= samples; ++i) { + double t = static_cast(i) / samples; + double ua = ea.u0 + t * (ea.u1 - ea.u0); + double va = ea.v0 + t * (ea.v1 - ea.v0); + double ub = eb.u0 + t * (eb.u1 - eb.u0); + double vb = eb.v0 + t * (eb.v1 - eb.v0); + + Point3D pa = surf_a.evaluate(ua, va); + Point3D pb = surf_b.evaluate(ub, vb); + double pos_err = (pa - pb).norm(); + report.position_errors.push_back(pos_err); + + ContinuityLevel level; + + if (pos_err > tol) { + level = ContinuityLevel::None; + } else { + Vector3D d1a = along_u_a ? surf_a.derivative_u(ua, va) : surf_a.derivative_v(ua, va); + Vector3D d1b = along_u_b ? surf_b.derivative_u(ub, vb) : surf_b.derivative_v(ub, vb); + + if (is_zero_vector(d1a) || is_zero_vector(d1b) || + !are_collinear_same_direction(d1a, d1b, tol)) { + level = ContinuityLevel::G0; + } else { + g1_count++; + + // G2 check via finite difference second derivative + Vector3D d2a, d2b; + double h = eps; + if (along_u_a) { + d2a = (surf_a.derivative_u(ua+h, va) - surf_a.derivative_u(ua-h, va)) / (2*h); + } else { + d2a = (surf_a.derivative_v(ua, va+h) - surf_a.derivative_v(ua, va-h)) / (2*h); + } + if (along_u_b) { + d2b = (surf_b.derivative_u(ub+h, vb) - surf_b.derivative_u(ub-h, vb)) / (2*h); + } else { + d2b = (surf_b.derivative_v(ub, vb+h) - surf_b.derivative_v(ub, vb-h)) / (2*h); + } + + if (!curvature_match(d1a, d2a, d1b, d2b, tol)) { + level = ContinuityLevel::G1; + } else { + g2_count++; + // G3 check + Vector3D d3a, d3b; + double h2 = eps * 2.0; + if (along_u_a) { + auto d2_fwd = (surf_a.derivative_u(ua+2*h2, va) - surf_a.derivative_u(ua,va)) / (2*h2); + auto d2_bwd = (surf_a.derivative_u(ua,va) - surf_a.derivative_u(ua-2*h2, va)) / (2*h2); + d3a = (d2_fwd - d2_bwd) / (2*h2); + } else { + auto d2_fwd = (surf_a.derivative_v(ua, va+2*h2) - surf_a.derivative_v(ua,va)) / (2*h2); + auto d2_bwd = (surf_a.derivative_v(ua,va) - surf_a.derivative_v(ua, va-2*h2)) / (2*h2); + d3a = (d2_fwd - d2_bwd) / (2*h2); + } + if (along_u_b) { + auto d2_fwd = (surf_b.derivative_u(ub+2*h2, vb) - surf_b.derivative_u(ub,vb)) / (2*h2); + auto d2_bwd = (surf_b.derivative_u(ub,vb) - surf_b.derivative_u(ub-2*h2, vb)) / (2*h2); + d3b = (d2_fwd - d2_bwd) / (2*h2); + } else { + auto d2_fwd = (surf_b.derivative_v(ub, vb+2*h2) - surf_b.derivative_v(ub,vb)) / (2*h2); + auto d2_bwd = (surf_b.derivative_v(ub,vb) - surf_b.derivative_v(ub, vb-2*h2)) / (2*h2); + d3b = (d2_fwd - d2_bwd) / (2*h2); + } + + if (!curvature_derivative_match(d1a, d2a, d3a, d1b, d2b, d3b, tol)) { + level = ContinuityLevel::G2; + } else { + level = ContinuityLevel::G3; + } + } + } + } + + report.per_sample.push_back(level); + if (static_cast(level) < static_cast(report.worst)) + report.worst = level; + if (static_cast(level) > static_cast(report.best)) + report.best = level; + } + + int total = static_cast(report.per_sample.size()); + if (total > 0) { + report.g1_ratio = static_cast(g1_count) / total; + report.g2_ratio = static_cast(g2_count) / total; + } + + return report; +} + +} // namespace vde::curves diff --git a/src/curves/surface_extension.cpp b/src/curves/surface_extension.cpp new file mode 100644 index 0000000..b7484df --- /dev/null +++ b/src/curves/surface_extension.cpp @@ -0,0 +1,636 @@ +#include "vde/curves/surface_extension.h" +#include "vde/curves/nurbs_operations.h" +#include "vde/curves/bspline_curve.h" +#include +#include +#include +#include +#include + +namespace vde::curves { + +using namespace vde::core; + +// =========================================================================== +// Helper: Greville abscissae for knot vector +// =========================================================================== +namespace { + +std::vector greville_abscissae(const std::vector& knots, int degree) { + int n = static_cast(knots.size()) - degree - 1; + std::vector greville(n); + for (int i = 0; i < n; ++i) { + double sum = 0.0; + for (int k = 1; k <= degree; ++k) + sum += knots[i + k]; + greville[i] = sum / degree; + } + return greville; +} + +/// Safe scalar-vector multiplication: returns Vector3D (not expression template) +inline Vector3D safe_mul(double scalar, const Vector3D& vec) { + return Vector3D(scalar * vec.x(), scalar * vec.y(), scalar * vec.z()); +} + +/// G1 linear extension implementation (replaces the call to nurbs_operations version) +NurbsSurface extend_surface_g1(const NurbsSurface& surf, int edge, double length) { + if (std::abs(length) < 1e-15) return surf; + if (length < 0.0) + throw std::invalid_argument("extend_surface: length must be non-negative"); + + const auto& cp = surf.control_points(); + const auto& w = surf.weights(); + auto ku = surf.knots_u(); + auto kv = surf.knots_v(); + int du = surf.degree_u(); + int dv = surf.degree_v(); + int nu = static_cast(cp.size()); + int nv = static_cast(cp.empty() ? 0 : cp[0].size()); + + if (nu == 0 || nv == 0) + throw std::invalid_argument("extend_surface: empty surface"); + + auto gu = greville_abscissae(ku, du); + auto gv = greville_abscissae(kv, dv); + + std::vector> new_cp; + std::vector> new_w; + std::vector new_ku, new_kv; + + switch (edge) { + case 0: { + double u_val = gu[0]; + std::vector new_row(nv); + std::vector new_row_w(nv, 1.0); + for (int j = 0; j < nv; ++j) { + Vector3D tan = surf.derivative_u(u_val, gv[j]); + double tlen = tan.norm(); + new_row[j] = Point3D( + cp[0][j].x() - (tlen > 1e-12 ? length * tan.x() / tlen : 0.0), + cp[0][j].y() - (tlen > 1e-12 ? length * tan.y() / tlen : 0.0), + cp[0][j].z() - (tlen > 1e-12 ? length * tan.z() / tlen : 0.0)); + if (!w.empty()) new_row_w[j] = w[0][j]; + } + new_cp.push_back(new_row); + for (int i = 0; i < nu; ++i) new_cp.push_back(cp[i]); + new_w = w; + if (!new_w.empty()) new_w.insert(new_w.begin(), new_row_w); + double dk = ku[du + nu - 1] - ku[du + nu - 2]; + new_ku = {ku[0] - dk}; + new_ku.insert(new_ku.end(), ku.begin(), ku.end()); + new_kv = kv; + break; + } + case 1: { + double u_val = gu[nu - 1]; + std::vector new_row(nv); + std::vector new_row_w(nv, 1.0); + for (int j = 0; j < nv; ++j) { + Vector3D tan = surf.derivative_u(u_val, gv[j]); + double tlen = tan.norm(); + new_row[j] = Point3D( + cp[nu - 1][j].x() + (tlen > 1e-12 ? length * tan.x() / tlen : 0.0), + cp[nu - 1][j].y() + (tlen > 1e-12 ? length * tan.y() / tlen : 0.0), + cp[nu - 1][j].z() + (tlen > 1e-12 ? length * tan.z() / tlen : 0.0)); + if (!w.empty()) new_row_w[j] = w[nu - 1][j]; + } + new_cp = cp; + new_cp.push_back(new_row); + new_w = w; + if (!new_w.empty()) new_w.push_back(new_row_w); + double dk = ku[du + nu - 1] - ku[du + nu - 2]; + new_ku = ku; + new_ku.push_back(ku.back() + dk); + new_kv = kv; + break; + } + case 2: { + double v_val = gv[0]; + for (int i = 0; i < nu; ++i) { + Vector3D tan = surf.derivative_v(gu[i], v_val); + double tlen = tan.norm(); + Point3D new_pt( + cp[i][0].x() - (tlen > 1e-12 ? length * tan.x() / tlen : 0.0), + cp[i][0].y() - (tlen > 1e-12 ? length * tan.y() / tlen : 0.0), + cp[i][0].z() - (tlen > 1e-12 ? length * tan.z() / tlen : 0.0)); + std::vector row = {new_pt}; + row.insert(row.end(), cp[i].begin(), cp[i].end()); + new_cp.push_back(row); + if (!w.empty()) { + std::vector row_w = {w[i][0]}; + row_w.insert(row_w.end(), w[i].begin(), w[i].end()); + new_w.push_back(row_w); + } + } + if (w.empty()) new_w = w; + double dk = kv[dv + nv - 1] - kv[dv + nv - 2]; + new_ku = ku; + new_kv = {kv[0] - dk}; + new_kv.insert(new_kv.end(), kv.begin(), kv.end()); + break; + } + case 3: { + double v_val = gv[nv - 1]; + for (int i = 0; i < nu; ++i) { + Vector3D tan = surf.derivative_v(gu[i], v_val); + double tlen = tan.norm(); + Point3D new_pt( + cp[i][nv - 1].x() + (tlen > 1e-12 ? length * tan.x() / tlen : 0.0), + cp[i][nv - 1].y() + (tlen > 1e-12 ? length * tan.y() / tlen : 0.0), + cp[i][nv - 1].z() + (tlen > 1e-12 ? length * tan.z() / tlen : 0.0)); + std::vector row = cp[i]; + row.push_back(new_pt); + new_cp.push_back(row); + if (!w.empty()) { + std::vector row_w = w[i]; + row_w.push_back(w[i][nv - 1]); + new_w.push_back(row_w); + } + } + if (w.empty()) new_w = w; + double dk = kv[dv + nv - 1] - kv[dv + nv - 2]; + new_ku = ku; + new_kv = kv; + new_kv.push_back(kv.back() + dk); + break; + } + default: + throw std::invalid_argument("extend_surface: edge must be 0-3"); + } + + return NurbsSurface(new_cp, new_ku, new_kv, new_w, du, dv); +} + +} // anonymous namespace + +// =========================================================================== +// extend_surface: G1 / G2 extension +// =========================================================================== +NurbsSurface extend_surface(const NurbsSurface& surface, + int edge_id, double distance, + Continuity continuity) { + if (std::abs(distance) < 1e-15) return surface; + if (distance < 0.0) + throw std::invalid_argument("extend_surface: distance must be non-negative"); + + const auto& cp = surface.control_points(); + const auto& w = surface.weights(); + auto ku = surface.knots_u(); + auto kv = surface.knots_v(); + int du = surface.degree_u(); + int dv = surface.degree_v(); + int nu = static_cast(cp.size()); + int nv = static_cast(cp.empty() ? 0 : cp[0].size()); + + if (nu == 0 || nv == 0) + throw std::invalid_argument("extend_surface: empty surface"); + + // For G1: use the local implementation + if (continuity == Continuity::G1) { + return extend_surface_g1(surface, edge_id, distance); + } + + // ── G2 continuation: add 2 rows/columns ── + auto gu = greville_abscissae(ku, du); + auto gv = greville_abscissae(kv, dv); + + std::vector> new_cp; + std::vector> new_w; + std::vector new_ku, new_kv; + + double dk = 0.0; // knot interval for extension + double half_dist = distance * 0.5; + + if (edge_id == 0 || edge_id == 1) { + // u-edge extension + dk = ku[du + nu - 1] - ku[du + nu - 2]; + double sign = (edge_id == 0) ? -1.0 : 1.0; + double u_val = (edge_id == 0) ? gu[0] : gu[nu - 1]; + int src_idx = (edge_id == 0) ? 0 : nu - 1; + + // --- Row 1: G1 tangent extrapolation --- + std::vector row1(nv); + std::vector row1_w(nv, 1.0); + for (int j = 0; j < nv; ++j) { + Vector3D tan = surface.derivative_u(u_val, gv[j]); + double tlen = tan.norm(); + if (tlen > 1e-12) tan = safe_mul(sign * half_dist / tlen, tan); + else tan = Vector3D::Zero(); + row1[j] = Point3D(cp[src_idx][j].x() + tan.x(), + cp[src_idx][j].y() + tan.y(), + cp[src_idx][j].z() + tan.z()); + if (!w.empty()) row1_w[j] = w[src_idx][j]; + } + + // --- Row 2: curvature-aware extrapolation --- + // Use the second derivative to adjust the 2nd row + std::vector row2(nv); + std::vector row2_w(nv, 1.0); + for (int j = 0; j < nv; ++j) { + Vector3D tan = surface.derivative_u(u_val, gv[j]); + double tlen = tan.norm(); + if (tlen > 1e-12) tan = safe_mul(sign * distance / tlen, tan); + else tan = Vector3D::Zero(); + + // Approximate curvature correction: extend further + // then adjust by curvature (2nd derivative info from adjacent CPs) + Point3D base_pt = Point3D( + cp[src_idx][j].x() + tan.x(), + cp[src_idx][j].y() + tan.y(), + cp[src_idx][j].z() + tan.z()); + + // Simple curvature estimate: use neighboring CP differences + if (nu >= 3 && du >= 2) { + int prev = (edge_id == 0) ? 1 : nu - 2; + Vector3D d1 = cp[src_idx][j] - cp[prev][j]; + Vector3D d2; + if (nu >= 4) { + int prev2 = (edge_id == 0) ? 2 : nu - 3; + d2 = cp[prev][j] - cp[prev2][j]; + } else { + d2 = d1; + } + // Curvature term = (d1 - d2) — second difference + Vector3D curv = d1 - d2; + base_pt = Point3D( + base_pt.x() + sign * curv.x() * 0.25, + base_pt.y() + sign * curv.y() * 0.25, + base_pt.z() + sign * curv.z() * 0.25); + } + row2[j] = base_pt; + if (!w.empty()) row2_w[j] = w[src_idx][j]; + } + + // Build new CP grid and knots + if (edge_id == 0) { + // prepend rows: row2, row1, then original CPs + new_cp.push_back(row2); + new_cp.push_back(row1); + for (int i = 0; i < nu; ++i) new_cp.push_back(cp[i]); + + if (!w.empty()) { + new_w.push_back(row2_w); + new_w.push_back(row1_w); + for (const auto& rw : w) new_w.push_back(rw); + } + + new_ku = {ku[0] - 2.0 * dk, ku[0] - dk}; + new_ku.insert(new_ku.end(), ku.begin(), ku.end()); + new_kv = kv; + } else { + // append rows: original CPs, then row1, row2 + new_cp = cp; + new_cp.push_back(row1); + new_cp.push_back(row2); + + new_w = w; + if (!w.empty()) { + new_w.push_back(row1_w); + new_w.push_back(row2_w); + } + + new_ku = ku; + new_ku.push_back(ku.back() + dk); + new_ku.push_back(ku.back() + 2.0 * dk); + new_kv = kv; + } + } else { + // v-edge extension (edge 2 or 3) + dk = kv[dv + nv - 1] - kv[dv + nv - 2]; + double sign = (edge_id == 2) ? -1.0 : 1.0; + double v_val = (edge_id == 2) ? gv[0] : gv[nv - 1]; + int src_idx = (edge_id == 2) ? 0 : nv - 1; + + // --- Column 1: G1 tangent extrapolation --- + // --- Column 2: curvature-aware extrapolation --- + for (int i = 0; i < nu; ++i) { + Vector3D tan = surface.derivative_v(gu[i], v_val); + double tlen = tan.norm(); + + Vector3D tan1 = Vector3D::Zero(); + Vector3D tan2 = Vector3D::Zero(); + if (tlen > 1e-12) { + tan1 = safe_mul(sign * half_dist / tlen, tan); + tan2 = safe_mul(sign * distance / tlen, tan); + } + + Point3D new_pt1(cp[i][src_idx].x() + tan1.x(), + cp[i][src_idx].y() + tan1.y(), + cp[i][src_idx].z() + tan1.z()); + Point3D new_pt2(cp[i][src_idx].x() + tan2.x(), + cp[i][src_idx].y() + tan2.y(), + cp[i][src_idx].z() + tan2.z()); + + // Curvature correction + if (nv >= 3 && dv >= 2) { + int prev = (edge_id == 2) ? 1 : nv - 2; + Vector3D d1 = cp[i][src_idx] - cp[i][prev]; + Vector3D d2; + if (nv >= 4) { + int prev2 = (edge_id == 2) ? 2 : nv - 3; + d2 = cp[i][prev] - cp[i][prev2]; + } else { + d2 = d1; + } + Vector3D curv = d1 - d2; + new_pt2 = Point3D( + new_pt2.x() + sign * curv.x() * 0.25, + new_pt2.y() + sign * curv.y() * 0.25, + new_pt2.z() + sign * curv.z() * 0.25); + } + + std::vector row; + std::vector row_w; + if (edge_id == 2) { + row = {new_pt2, new_pt1}; + row.insert(row.end(), cp[i].begin(), cp[i].end()); + if (!w.empty()) { + row_w = {w[i][src_idx], w[i][src_idx]}; + row_w.insert(row_w.end(), w[i].begin(), w[i].end()); + } + } else { + row = cp[i]; + row.push_back(new_pt1); + row.push_back(new_pt2); + if (!w.empty()) { + row_w = w[i]; + row_w.push_back(w[i][src_idx]); + row_w.push_back(w[i][src_idx]); + } + } + new_cp.push_back(row); + if (!w.empty()) new_w.push_back(row_w); + } + + new_ku = ku; + if (edge_id == 2) { + new_kv = {kv[0] - 2.0 * dk, kv[0] - dk}; + new_kv.insert(new_kv.end(), kv.begin(), kv.end()); + } else { + new_kv = kv; + new_kv.push_back(kv.back() + dk); + new_kv.push_back(kv.back() + 2.0 * dk); + } + } + + if (new_w.empty()) new_w = {}; + // Determine weights: if original had weights, use new_w; else empty + std::vector> final_w; + bool has_weights = !w.empty(); + if (has_weights && !new_w.empty()) final_w = new_w; + + return NurbsSurface(new_cp, new_ku, new_kv, final_w, + (edge_id == 0 || edge_id == 1) ? du : du, + (edge_id == 0 || edge_id == 1) ? dv : dv); +} + +// =========================================================================== +// n_sided_fill: N-sided hole filling with Coons patch generalization +// =========================================================================== +NurbsSurface n_sided_fill(const std::vector& boundary_curves, + Continuity continuity) { + size_t N = boundary_curves.size(); + if (N < 3) return NurbsSurface(); // degenerate + + // ── Step 1: Compute center point ── + Point3D center(0, 0, 0); + for (const auto& c : boundary_curves) { + center += c.evaluate(0.5); + } + center /= static_cast(N); + + // ── Step 2: Compute corner points (curve endpoints) ── + std::vector corners(N); + for (size_t i = 0; i < N; ++i) { + corners[i] = boundary_curves[i].evaluate(0.0); + } + + // ── Step 3: Sample boundary curves ── + int samples_per_curve = (continuity == Continuity::G2) ? 8 : 6; + int total_v = static_cast(N) * samples_per_curve; + + // Number of radial rings + int num_rings = (continuity == Continuity::G2) ? 4 : 3; + + std::vector> grid(num_rings, + std::vector(total_v)); + std::vector> weight_grid(num_rings, + std::vector(total_v, 1.0)); + + // ── Step 4: Fill rings ── + // Innermost ring → center + for (int j = 0; j < total_v; ++j) { + grid[0][j] = center; + } + + // Intermediate rings → radial interpolation from center to boundary + for (int ring = 1; ring < num_rings - 1; ++ring) { + double t = static_cast(ring) / (num_rings - 1); + for (size_t i = 0; i < N; ++i) { + for (int j = 0; j < samples_per_curve; ++j) { + double param = static_cast(j) / samples_per_curve; + Point3D bpt = boundary_curves[i].evaluate(param); + int idx = static_cast(i) * samples_per_curve + j; + + // Blend: linear for G1, smoothstep for G2 + double blend; + if (continuity == Continuity::G2) { + // Smoothstep: 3t² - 2t³ for smoother transition + double s = t; + blend = s * s * (3.0 - 2.0 * s); + } else { + blend = t; + } + grid[ring][idx] = Point3D( + center.x() + (1.0 - blend) * 0.0 + blend * (bpt.x() - center.x()), + center.y() + (1.0 - blend) * 0.0 + blend * (bpt.y() - center.y()), + center.z() + (1.0 - blend) * 0.0 + blend * (bpt.z() - center.z())); + // Simplified: grid[ring][idx] = center + blend * (bpt - center) + grid[ring][idx] = Point3D( + center.x() + blend * (bpt.x() - center.x()), + center.y() + blend * (bpt.y() - center.y()), + center.z() + blend * (bpt.z() - center.z())); + } + } + } + + // Outer ring → boundary points + for (size_t i = 0; i < N; ++i) { + for (int j = 0; j < samples_per_curve; ++j) { + double param = static_cast(j) / samples_per_curve; + int idx = static_cast(i) * samples_per_curve + j; + grid[num_rings - 1][idx] = boundary_curves[i].evaluate(param); + } + } + + // ── Step 5: Build knot vectors ── + // u-direction (radial): linear interpolation + std::vector ku(num_rings + 2); // nu + degree_u + 1 + for (int i = 0; i < static_cast(ku.size()); ++i) { + ku[i] = static_cast(i) / (ku.size() - 1); + } + + // v-direction (circumferential): B-spline with periodic-like clamping + // degree 2 in v, total_v control points + int deg_v = 2; + std::vector kv(total_v + deg_v + 1); + for (int i = 0; i < static_cast(kv.size()); ++i) { + kv[i] = static_cast(i) / (total_v + deg_v); + } + + // u-direction degree = 1 (linear in radial direction) + return NurbsSurface(grid, ku, kv, weight_grid, 1, deg_v); +} + +// =========================================================================== +// blend_surfaces: auto-detect common boundary +// =========================================================================== +NurbsSurface blend_surfaces(const NurbsSurface& surf_a, + const NurbsSurface& surf_b, + double radius) { + if (radius < 1e-15) return NurbsSurface(); + + const auto& cp_a = surf_a.control_points(); + const auto& cp_b = surf_b.control_points(); + int nu_a = static_cast(cp_a.size()); + int nv_a = static_cast(cp_a.empty() ? 0 : cp_a[0].size()); + int nu_b = static_cast(cp_b.size()); + int nv_b = static_cast(cp_b.empty() ? 0 : cp_b[0].size()); + + if (nu_a == 0 || nv_a == 0 || nu_b == 0 || nv_b == 0) + return NurbsSurface(); + + auto ku_a = surf_a.knots_u(), kv_a = surf_a.knots_v(); + auto ku_b = surf_b.knots_u(), kv_b = surf_b.knots_v(); + int du_a = surf_a.degree_u(), dv_a = surf_a.degree_v(); + int du_b = surf_b.degree_u(), dv_b = surf_b.degree_v(); + + auto gu_a = greville_abscissae(ku_a, du_a); + auto gv_a = greville_abscissae(kv_a, dv_a); + auto gu_b = greville_abscissae(ku_b, du_b); + auto gv_b = greville_abscissae(kv_b, dv_b); + + // Extract boundary curves for all 4 edges of each surface + // and find the common boundary. + struct EdgeInfo { + int edge_id; + std::vector endpoints; // 2 endpoints of the edge curve + NurbsCurve curve; + }; + + auto get_edge_info = [](const NurbsSurface& surf, int edge_id) -> EdgeInfo { + EdgeInfo info; + info.edge_id = edge_id; + info.curve = extract_boundary_curve(surf, edge_id); + double t0 = info.curve.domain().first; + double t1 = info.curve.domain().second; + info.endpoints = {info.curve.evaluate(t0), info.curve.evaluate(t1)}; + return info; + }; + + std::vector edges_a, edges_b; + for (int e = 0; e < 4; ++e) { + edges_a.push_back(get_edge_info(surf_a, e)); + edges_b.push_back(get_edge_info(surf_b, e)); + } + + // Find matching edge pair: endpoints must be close (within tolerance) + int best_edge_a = -1, best_edge_b = -1; + double best_dist = 1e30; + const double match_tol = 1e-6; + + for (int ea = 0; ea < 4; ++ea) { + for (int eb = 0; eb < 4; ++eb) { + double d1 = (edges_a[ea].endpoints[0] - edges_b[eb].endpoints[0]).norm(); + double d2 = (edges_a[ea].endpoints[1] - edges_b[eb].endpoints[1]).norm(); + double d3 = (edges_a[ea].endpoints[0] - edges_b[eb].endpoints[1]).norm(); + double d4 = (edges_a[ea].endpoints[1] - edges_b[eb].endpoints[0]).norm(); + + double d_forward = std::max(d1, d2); + double d_reverse = std::max(d3, d4); + double d = std::min(d_forward, d_reverse); + + if (d < match_tol && d < best_dist) { + best_dist = d; + best_edge_a = ea; + best_edge_b = eb; + } + } + } + + if (best_edge_a < 0 || best_edge_b < 0) { + // No common boundary found → return degenerate surface + return NurbsSurface(); + } + + // Use the existing blend_surfaces from nurbs_operations with explicit edge IDs + return ::vde::curves::blend_surfaces(surf_a, surf_b, + best_edge_a, best_edge_b, radius); +} + +// =========================================================================== +// is_g1_continuous +// =========================================================================== +bool is_g1_continuous(const NurbsSurface& surf_a, const NurbsSurface& surf_b, + int edge, int samples) { + const auto& ku_a = surf_a.knots_u(); + const auto& kv_a = surf_a.knots_v(); + const auto& ku_b = surf_b.knots_u(); + const auto& kv_b = surf_b.knots_v(); + + double u_min_a = ku_a[surf_a.degree_u()]; + double u_max_a = ku_a[ku_a.size() - surf_a.degree_u() - 1]; + double v_min_a = kv_a[surf_a.degree_v()]; + double v_max_a = kv_a[kv_a.size() - surf_a.degree_v() - 1]; + double u_min_b = ku_b[surf_b.degree_u()]; + double u_max_b = ku_b[ku_b.size() - surf_b.degree_u() - 1]; + double v_min_b = kv_b[surf_b.degree_v()]; + double v_max_b = kv_b[kv_b.size() - surf_b.degree_v() - 1]; + + const double tol_pos = 1e-6; + const double tol_norm = 1e-3; + + bool is_u_edge_a = (edge == 0 || edge == 1); + double u_fix_a = (edge == 0) ? u_min_a : ((edge == 1) ? u_max_a : 0.0); + double v_fix_a = (edge == 2) ? v_min_a : ((edge == 3) ? v_max_a : 0.0); + + int edge_b; + switch (edge) { + case 0: edge_b = 1; break; + case 1: edge_b = 0; break; + case 2: edge_b = 3; break; + case 3: edge_b = 2; break; + default: return false; + } + bool is_u_edge_b = (edge_b == 0 || edge_b == 1); + double u_fix_b = (edge_b == 0) ? u_min_b : ((edge_b == 1) ? u_max_b : 0.0); + double v_fix_b = (edge_b == 2) ? v_min_b : ((edge_b == 3) ? v_max_b : 0.0); + + for (int s = 0; s < samples; ++s) { + double t = static_cast(s) / std::max(1, samples - 1); + + double u_a, v_a; + if (is_u_edge_a) { u_a = u_fix_a; v_a = v_min_a + t * (v_max_a - v_min_a); } + else { u_a = u_min_a + t * (u_max_a - u_min_a); v_a = v_fix_a; } + + double u_b, v_b; + if (is_u_edge_b) { u_b = u_fix_b; v_b = v_min_b + t * (v_max_b - v_min_b); } + else { u_b = u_min_b + t * (u_max_b - u_min_b); v_b = v_fix_b; } + + // G0 check + Point3D pa = surf_a.evaluate(u_a, v_a); + Point3D pb = surf_b.evaluate(u_b, v_b); + if ((pa - pb).norm() > tol_pos) return false; + + // G1 check: tangent plane continuity (normals should align) + Vector3D Na = surf_a.normal(u_a, v_a); + Vector3D Nb = surf_b.normal(u_b, v_b); + double norm_dot = std::abs(Na.dot(Nb)); + if (std::abs(norm_dot - 1.0) > tol_norm) return false; + } + + return true; +} + +} // namespace vde::curves diff --git a/tests/brep/CMakeLists.txt b/tests/brep/CMakeLists.txt index b762af9..a10c686 100644 --- a/tests/brep/CMakeLists.txt +++ b/tests/brep/CMakeLists.txt @@ -31,3 +31,4 @@ add_vde_test(test_parallel_boolean) add_vde_test(test_feature_recognition) add_vde_test(test_v5_features) add_vde_test(test_v5_1) +add_vde_test(test_advanced_blend) diff --git a/tests/brep/test_advanced_blend.cpp b/tests/brep/test_advanced_blend.cpp new file mode 100644 index 0000000..af242ab --- /dev/null +++ b/tests/brep/test_advanced_blend.cpp @@ -0,0 +1,158 @@ +#include +#include "vde/brep/advanced_blend.h" +#include "vde/brep/modeling.h" +#include "vde/brep/brep_validate.h" +#include "vde/brep/brep.h" +#include + +using namespace vde::brep; +using namespace vde::core; + +// ═══════════════════════════════════════════════ +// variable_radius_blend 测试 (4项) +// ═══════════════════════════════════════════════ + +TEST(VariableRadiusBlendTest, BoxEdgeGrowRadius) { + auto box = make_box(4, 4, 4); + auto result = variable_radius_blend(box, 0, 0.3, 0.8, 16); + // Blend adds blend-strip faces → more faces than original + EXPECT_GT(result.num_faces(), box.num_faces()); + EXPECT_GT(result.num_edges(), box.num_edges()); + EXPECT_GT(result.num_vertices(), 0u); + EXPECT_TRUE(result.is_valid()); +} + +TEST(VariableRadiusBlendTest, BoxEdgeShrinkRadius) { + auto box = make_box(4, 4, 4); + auto result = variable_radius_blend(box, 0, 0.8, 0.3, 12); + EXPECT_GT(result.num_faces(), box.num_faces()); + EXPECT_TRUE(result.is_valid()); +} + +TEST(VariableRadiusBlendTest, BoxEdgeSymmetric) { + auto box = make_box(4, 4, 4); + // Equal start and end radius → constant-radius blend + auto r1 = variable_radius_blend(box, 0, 0.5, 0.5, 8); + auto r2 = variable_radius_blend(box, 0, 0.5, 0.5, 16); + EXPECT_GT(r1.num_faces(), box.num_faces()); + // More samples = more blend strips = more faces + EXPECT_GT(r2.num_faces(), r1.num_faces()); + EXPECT_TRUE(r1.is_valid()); + EXPECT_TRUE(r2.is_valid()); +} + +TEST(VariableRadiusBlendTest, SmallRadiusNearZero) { + auto box = make_box(4, 4, 4); + auto result = variable_radius_blend(box, 1, 0.05, 0.5, 8); + EXPECT_GT(result.num_faces(), box.num_faces()); + EXPECT_TRUE(result.is_valid()); +} + +// ═══════════════════════════════════════════════ +// multi_face_blend 测试 (3项) +// ═══════════════════════════════════════════════ + +TEST(MultiFaceBlendTest, TwoFaces) { + auto box = make_box(4, 4, 4); + auto result = multi_face_blend(box, {0, 1}, 0.5); + EXPECT_GT(result.num_vertices(), 0u); + EXPECT_TRUE(result.is_valid()); +} + +TEST(MultiFaceBlendTest, ThreeFacesCorner) { + // Blend 3 faces meeting at a corner — spherical corner fill added + auto box = make_box(4, 4, 4); + auto result = multi_face_blend(box, {0, 2, 4}, 0.4); + EXPECT_GT(result.num_vertices(), 0u); + EXPECT_TRUE(result.is_valid()); +} + +TEST(MultiFaceBlendTest, AllFaces) { + auto box = make_box(4, 4, 4); + std::vector all_faces; + for (size_t i = 0; i < box.num_faces(); ++i) + all_faces.push_back(static_cast(i)); + auto result = multi_face_blend(box, all_faces, 0.3); + EXPECT_GT(result.num_vertices(), 0u); + EXPECT_TRUE(result.is_valid()); +} + +// ═══════════════════════════════════════════════ +// rolling_ball_blend 测试 (4项) +// ═══════════════════════════════════════════════ + +TEST(RollingBallBlendTest, StraightEdge) { + auto box = make_box(4, 4, 4); + auto result = rolling_ball_blend(box, 0, 0.5, 20); + EXPECT_GT(result.num_faces(), box.num_faces()); + EXPECT_TRUE(result.is_valid()); +} + +TEST(RollingBallBlendTest, ConvexCorner) { + // Box edge is a convex corner + auto box = make_box(4, 4, 4); + auto result = rolling_ball_blend(box, 3, 0.6, 16); + EXPECT_GT(result.num_faces(), box.num_faces()); + EXPECT_TRUE(result.is_valid()); +} + +TEST(RollingBallBlendTest, DifferentEdges) { + auto box = make_box(4, 4, 4); + // Test rolling ball on first 4 edges + for (int ei = 0; ei < 4; ++ei) { + auto result = rolling_ball_blend(box, ei, 0.3, 10); + EXPECT_GT(result.num_faces(), box.num_faces()) << "edge " << ei; + EXPECT_TRUE(result.is_valid()) << "edge " << ei; + } +} + +TEST(RollingBallBlendTest, VaryingSamples) { + auto box = make_box(4, 4, 4); + auto r1 = rolling_ball_blend(box, 0, 0.5, 8); + auto r2 = rolling_ball_blend(box, 0, 0.5, 20); + // More samples => more blend strip faces + EXPECT_GT(r2.num_faces(), r1.num_faces()); + EXPECT_TRUE(r1.is_valid()); + EXPECT_TRUE(r2.is_valid()); +} + +// ═══════════════════════════════════════════════ +// face_face_blend 测试 (4项) +// ═══════════════════════════════════════════════ + +TEST(FaceFaceBlendTest, AdjacentFaces) { + auto box = make_box(4, 4, 4); + // face 0 (-Z) and face 2 (-Y) share an edge + auto result = face_face_blend(box, 0, 2, 0.5, 16); + EXPECT_GT(result.num_faces(), box.num_faces()); + EXPECT_TRUE(result.is_valid()); +} + +TEST(FaceFaceBlendTest, DifferentAdjacentPair) { + auto box = make_box(4, 4, 4); + auto result = face_face_blend(box, 1, 4, 0.4, 12); + EXPECT_GT(result.num_faces(), box.num_faces()); + EXPECT_TRUE(result.is_valid()); +} + +TEST(FaceFaceBlendTest, MultiplePairs) { + auto box = make_box(4, 4, 4); + auto r1 = face_face_blend(box, 0, 2, 0.3, 8); + auto r2 = face_face_blend(box, 1, 4, 0.4, 10); + auto r3 = face_face_blend(box, 3, 5, 0.35, 12); + EXPECT_GT(r1.num_faces(), box.num_faces()); + EXPECT_GT(r2.num_faces(), box.num_faces()); + EXPECT_GT(r3.num_faces(), box.num_faces()); + EXPECT_TRUE(r1.is_valid()); + EXPECT_TRUE(r2.is_valid()); + EXPECT_TRUE(r3.is_valid()); +} + +TEST(FaceFaceBlendTest, NonSharedEdgeFallback) { + // Opposite faces — no shared edge, falls back gracefully + auto box = make_box(4, 4, 4); + auto result = face_face_blend(box, 0, 1, 0.2, 8); + // Should return a valid model (may be original if no common edge) + EXPECT_GT(result.num_faces(), 0u); + EXPECT_TRUE(result.is_valid()); +} diff --git a/tests/curves/CMakeLists.txt b/tests/curves/CMakeLists.txt index d9627f8..c99c722 100644 --- a/tests/curves/CMakeLists.txt +++ b/tests/curves/CMakeLists.txt @@ -3,3 +3,6 @@ add_vde_test(test_nurbs) add_vde_test(test_nurbs_operations) add_vde_test(test_surface_intersection) add_vde_test(test_parallel_intersection) +add_vde_test(test_surface_continuity) +add_vde_test(test_surface_analysis) +add_vde_test(test_surface_extension) diff --git a/tests/curves/test_surface_analysis.cpp b/tests/curves/test_surface_analysis.cpp new file mode 100644 index 0000000..f7021e7 --- /dev/null +++ b/tests/curves/test_surface_analysis.cpp @@ -0,0 +1,181 @@ +#include +#include "vde/curves/surface_analysis.h" +#include "vde/curves/nurbs_curve.h" +#include "vde/curves/nurbs_surface.h" +#include "vde/brep/modeling.h" +#include + +using namespace vde::curves; +using namespace vde::core; +using namespace vde::brep; + +// --------------------------------------------------------------------------- +// 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 simple NURBS curve (line) +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 (parabola-like) +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: CurvatureMap +// --------------------------------------------------------------------------- + +TEST(SurfaceAnalysisTest, CurvatureMap_PlaneIsDevelopable) { + auto surf = plane_surface(); + auto cm = curvature_map(surf, 20, 20); + + EXPECT_EQ(21, static_cast(cm.grid.size())); + EXPECT_EQ(21, static_cast(cm.grid[0].size())); + + // Plane: Gaussian curvature ≈ 0, Mean curvature ≈ 0 + EXPECT_NEAR(cm.min_gaussian, 0.0, 1e-9); + EXPECT_NEAR(cm.max_gaussian, 0.0, 1e-9); + EXPECT_NEAR(cm.min_mean, 0.0, 1e-9); + EXPECT_NEAR(cm.max_mean, 0.0, 1e-9); + + EXPECT_TRUE(cm.is_developable()); +} + +TEST(SurfaceAnalysisTest, CurvatureMap_GridDimensions) { + auto surf = plane_surface(); + auto cm = curvature_map(surf, 5, 8); + + EXPECT_EQ(6, static_cast(cm.grid.size())); // res_u+1 + EXPECT_EQ(9, static_cast(cm.grid[0].size())); // res_v+1 + EXPECT_EQ(5, cm.res_u); + EXPECT_EQ(8, cm.res_v); +} + +// --------------------------------------------------------------------------- +// Test: Zebra Stripe +// --------------------------------------------------------------------------- + +TEST(SurfaceAnalysisTest, ZebraStripe_OutputDimensions) { + auto surf = plane_surface(); + auto zs = zebra_stripe(surf, Vector3D(1, 0, 0), 15, 10); + + EXPECT_EQ(16, static_cast(zs.intensity.size())); + EXPECT_EQ(11, static_cast(zs.intensity[0].size())); + EXPECT_EQ(15, zs.res_u); + EXPECT_EQ(10, zs.res_v); + + // Check light direction is normalized + EXPECT_NEAR(zs.light_direction.norm(), 1.0, 1e-12); +} + +TEST(SurfaceAnalysisTest, ZebraStripe_IntensityRange) { + auto surf = plane_surface(); + auto zs = zebra_stripe(surf, Vector3D(0, 0, 1), 10, 10); + + for (size_t i = 0; i < zs.intensity.size(); ++i) { + for (size_t j = 0; j < zs.intensity[i].size(); ++j) { + double val = zs.intensity[i][j]; + EXPECT_GE(val, 0.0); + EXPECT_LE(val, 1.0); + } + } +} + +// --------------------------------------------------------------------------- +// Test: Curvature Comb +// --------------------------------------------------------------------------- + +TEST(SurfaceAnalysisTest, CurvatureComb_StraightLine) { + auto c = line(Point3D(0, 0, 0), Point3D(10, 0, 0)); + auto comb = curvature_comb(c, 20); + + EXPECT_EQ(21, static_cast(comb.points.size())); + // Straight line has zero curvature + EXPECT_NEAR(comb.min_curvature, 0.0, 1e-12); + EXPECT_NEAR(comb.max_curvature, 0.0, 1e-12); +} + +TEST(SurfaceAnalysisTest, CurvatureComb_NonzeroCurvature) { + auto c = cubic_curve(); + auto comb = curvature_comb(c, 100); + + EXPECT_EQ(101, static_cast(comb.points.size())); + // Cubic curve has non-zero curvature + EXPECT_GT(comb.max_curvature, 0.0); + // Scale should be auto-computed + EXPECT_GT(comb.scale, 0.0); +} + +// --------------------------------------------------------------------------- +// Test: Deviation Analysis +// --------------------------------------------------------------------------- + +TEST(SurfaceAnalysisTest, DeviationAnalysis_IdenticalPlanes) { + auto s1 = plane_surface(); + auto s2 = plane_surface(); + auto result = deviation_analysis(s1, s2, 10); + + EXPECT_EQ(100u, result.sample_count); + // Identical planes: deviation ≈ 0 + EXPECT_NEAR(result.rms, 0.0, 1e-6); + EXPECT_NEAR(result.mean_absolute, 0.0, 1e-6); +} + +TEST(SurfaceAnalysisTest, DeviationAnalysis_OffsetPlanes) { + auto s1 = plane_surface(); + + // Same plane shifted up by 1 in Z + std::vector> grid2 = { + {Point3D(0,0,1), Point3D(0,1,1)}, + {Point3D(1,0,1), Point3D(1,1,1)} + }; + auto s2 = NurbsSurface(grid2, {0,0,1,1}, {0,0,1,1}, + {{1,1},{1,1}}, 1, 1); + + auto result = deviation_analysis(s1, s2, 10); + + EXPECT_NEAR(result.max_positive, 1.0, 0.05); + EXPECT_NEAR(result.rms, 1.0, 0.05); +} + +// --------------------------------------------------------------------------- +// Test: Draft Face Angle +// --------------------------------------------------------------------------- + +TEST(SurfaceAnalysisTest, DraftAngle_HorizontalPlaneVerticalPull) { + auto surf = plane_surface(); + // Horizontal plane (normal = +Z), vertical pull (+Z) → angle = -π/2 + // (face normal parallel to pull: face is perpendicular to opening) + double angle = draft_face_angle(surf, 0.5, 0.5, Vector3D(0, 0, 1)); + EXPECT_NEAR(angle, -M_PI_2, 1e-6); +} + +TEST(SurfaceAnalysisTest, DraftAngle_PullPerpendicularToNormal) { + auto surf = plane_surface(); + // Horizontal plane (normal = +Z), pull in X → face parallel to pull → angle ≈ 0 + double angle = draft_face_angle(surf, 0.5, 0.5, Vector3D(1, 0, 0)); + EXPECT_NEAR(angle, 0.0, 1e-6); +} + +TEST(SurfaceAnalysisTest, DraftAngle_BrepModel) { + auto box = make_box(10, 10, 10); + // Top face has normal ≈ +Z, pull +Z → angle ≈ -π/2 + double angle = draft_face_angle(box, 0, Vector3D(0, 0, 1), 5); + EXPECT_NEAR(std::abs(angle), M_PI_2, 0.01); +} diff --git a/tests/curves/test_surface_continuity.cpp b/tests/curves/test_surface_continuity.cpp new file mode 100644 index 0000000..70a9818 --- /dev/null +++ b/tests/curves/test_surface_continuity.cpp @@ -0,0 +1,195 @@ +#include +#include "vde/curves/surface_continuity.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 degree-1 NURBS line from a to b +static NurbsCurve line(const Point3D& a, const Point3D& b) { + return NurbsCurve({a, b}, {0, 0, 1, 1}, {1, 1}, 1); +} + +/// 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 simple NURBS surface from a 2×2 grid +static NurbsSurface simple_surface(const std::vector>& grid, + int du=1, int dv=1) { + int rows = static_cast(grid.size()); + int cols = static_cast(grid.empty() ? 0 : grid[0].size()); + std::vector ku(rows + du + 1, 0.0); + std::vector kv(cols + dv + 1, 0.0); + for (int i = 0; i <= du; ++i) ku[i] = 0.0; + for (size_t i = du + 1; i < ku.size() - du - 1; ++i) + ku[i] = static_cast(i - du) / (ku.size() - 2*du - 1); + for (size_t i = ku.size() - du - 1; i < ku.size(); ++i) ku[i] = 1.0; + for (int i = 0; i <= dv; ++i) kv[i] = 0.0; + for (size_t i = dv + 1; i < kv.size() - dv - 1; ++i) + kv[i] = static_cast(i - dv) / (kv.size() - 2*dv - 1); + for (size_t i = kv.size() - dv - 1; i < kv.size(); ++i) kv[i] = 1.0; + + std::vector> w(rows, std::vector(cols, 1.0)); + return NurbsSurface(grid, ku, kv, w, du, dv); +} + +// --------------------------------------------------------------------------- +// Test: to_string +// --------------------------------------------------------------------------- + +TEST(ContinuityTest, ToString) { + EXPECT_STREQ("G0", to_string(ContinuityLevel::G0)); + EXPECT_STREQ("G1", to_string(ContinuityLevel::G1)); + EXPECT_STREQ("G2", to_string(ContinuityLevel::G2)); + EXPECT_STREQ("G3", to_string(ContinuityLevel::G3)); + EXPECT_STREQ("None", to_string(ContinuityLevel::None)); +} + +// --------------------------------------------------------------------------- +// Test: G0 — Position continuity on lines +// --------------------------------------------------------------------------- + +TEST(ContinuityTest, G0_SameCurveContinuous) { + // Two identical line segments should be G1 (same tangent direction) + auto a = line(Point3D(0, 0, 0), Point3D(1, 0, 0)); + auto b = line(Point3D(1, 0, 0), Point3D(2, 0, 0)); + auto level = continuity_type(a, b); + EXPECT_GE(static_cast(level), static_cast(ContinuityLevel::G1)); +} + +TEST(ContinuityTest, G0_EndpointsMatch) { + auto a = line(Point3D(0, 0, 0), Point3D(1, 0, 0)); + auto b = line(Point3D(1, 0, 0), Point3D(1, 1, 0)); + auto level = continuity_type(a, b); + EXPECT_GE(static_cast(level), static_cast(ContinuityLevel::G0)); +} + +TEST(ContinuityTest, G0_NoConnection) { + auto a = line(Point3D(0, 0, 0), Point3D(1, 0, 0)); + auto b = line(Point3D(5, 5, 5), Point3D(6, 6, 6)); + auto level = continuity_type(a, b); + EXPECT_EQ(ContinuityLevel::None, level); +} + +// --------------------------------------------------------------------------- +// Test: G1 — Tangent continuity +// --------------------------------------------------------------------------- + +TEST(ContinuityTest, G1_CollinearLines) { + // Two collinear lines → G0 position, G1 tangent + auto a = line(Point3D(0, 0, 0), Point3D(1, 0, 0)); + auto b = line(Point3D(1, 0, 0), Point3D(2, 0, 0)); + // Check at the connection point with explicit t values + auto level = continuity_type(a, 1.0, b, 0.0); + EXPECT_GE(static_cast(level), static_cast(ContinuityLevel::G1)); +} + +TEST(ContinuityTest, G1_SharpCornerG0Only) { + // Two lines meeting at a right angle → only G0 + auto a = line(Point3D(0, 0, 0), Point3D(1, 0, 0)); + auto b = line(Point3D(1, 0, 0), Point3D(1, 1, 0)); + auto level = continuity_type(a, 1.0, b, 0.0); + EXPECT_EQ(ContinuityLevel::G0, level); +} + +// --------------------------------------------------------------------------- +// Test: G2 — Curvature continuity on curves +// --------------------------------------------------------------------------- + +TEST(ContinuityTest, G2_StraightLinesG3) { + // Two collinear straight lines have zero curvature everywhere → G3 + auto a = line(Point3D(0, 0, 0), Point3D(1, 0, 0)); + auto b = line(Point3D(1, 0, 0), Point3D(3, 0, 0)); + auto level = continuity_type(a, 1.0, b, 0.0); + // Straight lines have d2=0, curvature=0, so they trivially match → G3 + EXPECT_EQ(ContinuityLevel::G3, level); +} + +TEST(ContinuityTest, G0_DisconnectedReturnsNone) { + auto a = line(Point3D(0, 0, 0), Point3D(1, 0, 0)); + auto b = line(Point3D(10, 10, 10), Point3D(11, 11, 11)); + auto level = continuity_type(a, 0.5, b, 0.5); + EXPECT_EQ(ContinuityLevel::None, level); +} + +// --------------------------------------------------------------------------- +// Test: Surface Continuity +// --------------------------------------------------------------------------- + +TEST(SurfaceContinuityTest, IdenticalPlanes) { + // Two identical flat planes → should be G3 + auto s1 = plane_surface(); + auto s2 = plane_surface(); + auto level = surface_continuity(s1, s2, 10); + EXPECT_GE(static_cast(level), static_cast(ContinuityLevel::G3)); +} + +TEST(SurfaceContinuityTest, CoplanarSurfacesEdgeMatch) { + // Two adjacent coplanar surfaces sharing an edge + std::vector> g1 = { + {Point3D(0,0,0), Point3D(0,1,0)}, + {Point3D(1,0,0), Point3D(1,1,0)} + }; + std::vector> g2 = { + {Point3D(1,0,0), Point3D(1,1,0)}, + {Point3D(2,0,0), Point3D(2,1,0)} + }; + auto s1 = simple_surface(g1); + auto s2 = simple_surface(g2); + auto level = surface_continuity(s1, s2, 10); + EXPECT_GE(static_cast(level), static_cast(ContinuityLevel::G3)); +} + +TEST(SurfaceContinuityTest, DisconnectedSurfaces) { + std::vector> g1 = { + {Point3D(0,0,0), Point3D(0,1,0)}, + {Point3D(1,0,0), Point3D(1,1,0)} + }; + std::vector> g2 = { + {Point3D(10,0,0), Point3D(10,1,0)}, + {Point3D(11,0,0), Point3D(11,1,0)} + }; + auto s1 = simple_surface(g1); + auto s2 = simple_surface(g2); + auto level = surface_continuity(s1, s2, 10); + EXPECT_EQ(ContinuityLevel::None, level); +} + +TEST(SurfaceContinuityTest, ContinuityReport) { + auto s1 = plane_surface(); + auto s2 = plane_surface(); + auto report = surface_continuity_report(s1, s2, 20); + EXPECT_GE(static_cast(report.worst), static_cast(ContinuityLevel::G3)); + EXPECT_NEAR(report.g1_ratio, 1.0, 0.01); + EXPECT_NEAR(report.g2_ratio, 1.0, 0.01); +} + +// --------------------------------------------------------------------------- +// Test: Edge parameter overload +// --------------------------------------------------------------------------- + +TEST(SurfaceContinuityTest, ExplicitEdgeParametersPlane) { + auto s1 = plane_surface(); + auto s2 = plane_surface(); + + // Same physical edge (u=1) on both identical surfaces + auto level = surface_continuity( + s1, s2, + {{1.0, 0.0}, {1.0, 1.0}}, // s1: u=1 edge, v 0→1 + {{1.0, 0.0}, {1.0, 1.0}}, // s2: same u=1 edge + 20); + EXPECT_GE(static_cast(level), static_cast(ContinuityLevel::G3)); +} diff --git a/tests/curves/test_surface_extension.cpp b/tests/curves/test_surface_extension.cpp new file mode 100644 index 0000000..a559e3d --- /dev/null +++ b/tests/curves/test_surface_extension.cpp @@ -0,0 +1,303 @@ +#include +#include "vde/curves/surface_extension.h" +#include "vde/curves/nurbs_operations.h" +#include "vde/curves/nurbs_curve.h" +#include "vde/curves/nurbs_surface.h" +#include + +using namespace vde::curves; +using namespace vde::core; + +// =========================================================================== +// Helpers +// =========================================================================== + +static NurbsCurve make_line_curve(const Point3D& a, const Point3D& b) { + return NurbsCurve({a, b}, {0, 0, 1, 1}, {1, 1}, 1); +} + +static NurbsSurface make_planar_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); +} + +static NurbsSurface make_degree2_surface() { + std::vector> grid = { + {Point3D(0, 0, 0), Point3D(0, 1, 0.3), Point3D(0, 2, 0)}, + {Point3D(1, 0, 0.3), Point3D(1, 1, 0.8), Point3D(1, 2, 0.3)}, + {Point3D(2, 0, 0), Point3D(2, 1, 0.3), Point3D(2, 2, 0)} + }; + return NurbsSurface(grid, + {0, 0, 0, 1, 1, 1}, {0, 0, 0, 1, 1, 1}, {}, 2, 2); +} + +// =========================================================================== +// extend_surface — G1 延伸 +// =========================================================================== + +TEST(SurfaceExtensionTest, ExtendG1_Umax_AddsRow) { + auto plane = make_planar_surface(); + auto extended = extend_surface(plane, 1, 0.5, Continuity::G1); + + // Should have one extra row of control points + auto [nu, nv] = extended.num_control_points(); + EXPECT_EQ(nu, 3); // original 2 + 1 new row + EXPECT_EQ(nv, 2); + + // Extended surface should evaluate to finite values throughout + for (double u = 0.0; u <= 1.0; u += 0.25) { + for (double v = 0.0; v <= 1.0; v += 0.25) { + Point3D p = extended.evaluate(u, v); + EXPECT_TRUE(std::isfinite(p.x())); + EXPECT_TRUE(std::isfinite(p.y())); + EXPECT_TRUE(std::isfinite(p.z())); + } + } +} + +TEST(SurfaceExtensionTest, ExtendG1_Vmin_AddsColumn) { + auto plane = make_planar_surface(); + auto extended = extend_surface(plane, 2, 1.0, Continuity::G1); + + auto [nu, nv] = extended.num_control_points(); + EXPECT_EQ(nu, 2); + EXPECT_EQ(nv, 3); // original 2 + 1 new column + + // All points should be finite + Point3D p = extended.evaluate(0.5, 0.5); + EXPECT_TRUE(std::isfinite(p.x())); + EXPECT_TRUE(std::isfinite(p.y())); + EXPECT_TRUE(std::isfinite(p.z())); +} + +TEST(SurfaceExtensionTest, ExtendG1_ZeroDistance) { + auto plane = make_planar_surface(); + auto same = extend_surface(plane, 0, 0.0, Continuity::G1); + + // Should return the same surface + Point3D p1 = plane.evaluate(0.5, 0.5); + Point3D p2 = same.evaluate(0.5, 0.5); + EXPECT_NEAR((p1 - p2).norm(), 0.0, 1e-10); +} + +// =========================================================================== +// extend_surface — G2 延伸 +// =========================================================================== + +TEST(SurfaceExtensionTest, ExtendG2_Umax_AddsTwoRows) { + auto plane = make_planar_surface(); + auto extended = extend_surface(plane, 1, 0.5, Continuity::G2); + + auto [nu, nv] = extended.num_control_points(); + EXPECT_EQ(nu, 4); // original 2 + 2 new rows + EXPECT_EQ(nv, 2); +} + +TEST(SurfaceExtensionTest, ExtendG2_PlanarStaysPlanar) { + auto plane = make_planar_surface(); + auto extended = extend_surface(plane, 1, 1.0, Continuity::G2); + + // All points should have z=0 (planar extension of XY plane) + for (double u = 0.0; u <= 1.0; u += 0.25) { + for (double v = 0.0; v <= 1.0; v += 0.25) { + Point3D p = extended.evaluate(u, v); + EXPECT_NEAR(p.z(), 0.0, 1e-6); + } + } +} + +TEST(SurfaceExtensionTest, ExtendG2_VmaxDegree2_AddsTwoCols) { + auto surf = make_degree2_surface(); + auto extended = extend_surface(surf, 3, 0.5, Continuity::G2); + + // Should have original 3 + 2 new columns = 5 in v + auto [nu, nv] = extended.num_control_points(); + EXPECT_EQ(nu, 3); + EXPECT_EQ(nv, 5); + + // Evaluate at valid parameters + Point3D p = extended.evaluate(0.5, 0.5); + EXPECT_TRUE(std::isfinite(p.x())); + EXPECT_TRUE(std::isfinite(p.y())); + EXPECT_TRUE(std::isfinite(p.z())); +} + +TEST(SurfaceExtensionTest, ExtendG1_PlanarPreservesFlatness) { + auto plane = make_planar_surface(); + auto extended = extend_surface(plane, 1, 0.5, Continuity::G1); + + // Extended surface should remain planar + for (double u = 0.0; u <= 1.0; u += 0.25) { + for (double v = 0.0; v <= 1.0; v += 0.25) { + Point3D p = extended.evaluate(u, v); + EXPECT_NEAR(p.z(), 0.0, 1e-6); + } + } +} + +// =========================================================================== +// n_sided_fill — N 边填充 +// =========================================================================== + +TEST(SurfaceExtensionTest, FillNSided_TriangleG2) { + auto l1 = make_line_curve(Point3D(0, 0, 0), Point3D(2, 0, 0)); + auto l2 = make_line_curve(Point3D(2, 0, 0), Point3D(1, 2, 0)); + auto l3 = make_line_curve(Point3D(1, 2, 0), Point3D(0, 0, 0)); + std::vector boundaries = {l1, l2, l3}; + + auto fill = n_sided_fill(boundaries, Continuity::G2); + + // Should produce a valid surface + Point3D mid = fill.evaluate(0.5, 0.5); + EXPECT_TRUE(std::isfinite(mid.x())); + EXPECT_TRUE(std::isfinite(mid.y())); + EXPECT_TRUE(std::isfinite(mid.z())); + + // Center should be roughly inside the triangle + EXPECT_GT(mid.x(), 0.0); + EXPECT_LT(mid.x(), 2.0); + EXPECT_GT(mid.y(), 0.0); + EXPECT_LT(mid.y(), 2.0); +} + +TEST(SurfaceExtensionTest, FillNSided_Quadrilateral) { + auto l1 = make_line_curve(Point3D(0, 0, 0), Point3D(3, 0, 0)); + auto l2 = make_line_curve(Point3D(3, 0, 0), Point3D(3, 2, 0)); + auto l3 = make_line_curve(Point3D(3, 2, 0), Point3D(0, 2, 0)); + auto l4 = make_line_curve(Point3D(0, 2, 0), Point3D(0, 0, 0)); + std::vector boundaries = {l1, l2, l3, l4}; + + auto fill = n_sided_fill(boundaries, Continuity::G1); + + // Center should be inside the rectangle + Point3D mid = fill.evaluate(0.5, 0.5); + EXPECT_GT(mid.x(), 0.5); + EXPECT_LT(mid.x(), 2.5); + EXPECT_GT(mid.y(), 0.0); + EXPECT_LT(mid.y(), 2.0); + EXPECT_NEAR(mid.z(), 0.0, 1e-6); +} + +TEST(SurfaceExtensionTest, FillNSided_Pentagon) { + // Regular pentagon in XY plane + std::vector boundaries; + std::vector pts; + for (int i = 0; i < 5; ++i) { + double angle = 2.0 * M_PI * i / 5.0 - M_PI / 2.0; + pts.push_back(Point3D(std::cos(angle), std::sin(angle), 0)); + } + for (int i = 0; i < 5; ++i) { + boundaries.push_back(make_line_curve(pts[i], pts[(i + 1) % 5])); + } + + auto fill = n_sided_fill(boundaries, Continuity::G2); + + // Center should be near origin + Point3D mid = fill.evaluate(0.5, 0.5); + EXPECT_NEAR(mid.x(), 0.0, 0.5); + EXPECT_NEAR(mid.y(), 0.0, 0.5); + EXPECT_NEAR(mid.z(), 0.0, 1e-6); +} + +TEST(SurfaceExtensionTest, FillNSided_TriangleG1) { + auto l1 = make_line_curve(Point3D(0, 0, 0), Point3D(2, 0, 0)); + auto l2 = make_line_curve(Point3D(2, 0, 0), Point3D(0, 2, 0)); + auto l3 = make_line_curve(Point3D(0, 2, 0), Point3D(0, 0, 0)); + std::vector boundaries = {l1, l2, l3}; + + auto fill = n_sided_fill(boundaries, Continuity::G1); + + // Fill surface should be finite and pass through boundary + for (double u = 0.0; u <= 1.0; u += 0.25) { + for (double v = 0.0; v <= 1.0; v += 0.25) { + Point3D p = fill.evaluate(u, v); + EXPECT_TRUE(std::isfinite(p.x())); + EXPECT_TRUE(std::isfinite(p.y())); + EXPECT_TRUE(std::isfinite(p.z())); + } + } +} + +TEST(SurfaceExtensionTest, FillNSided_DegenerateEdge) { + // Only 2 curves → should return degenerate + auto l1 = make_line_curve(Point3D(0, 0, 0), Point3D(1, 0, 0)); + auto l2 = make_line_curve(Point3D(1, 0, 0), Point3D(0, 1, 0)); + std::vector boundaries = {l1, l2}; + + auto fill = n_sided_fill(boundaries, Continuity::G1); + auto [nu, nv] = fill.num_control_points(); + // With < 3 curves, returns empty surface + EXPECT_TRUE(nu == 0 || nu == 1); +} + +// =========================================================================== +// blend_surfaces — 两面过渡（自动检测公共边） +// =========================================================================== + +TEST(SurfaceExtensionTest, BlendSurfacesAutoDetect_Adjacent) { + auto plane1 = make_planar_surface(); + std::vector> grid2 = { + {Point3D(1, 0, 0), Point3D(1, 1, 0)}, + {Point3D(2, 0, 0), Point3D(2, 1, 0)} + }; + NurbsSurface plane2(grid2, {0, 0, 1, 1}, {0, 0, 1, 1}, {}, 1, 1); + + auto blend = blend_surfaces(plane1, plane2, 0.3); + + // Should produce a valid blend surface + auto [nu, nv] = blend.num_control_points(); + EXPECT_GE(nu, 2); + EXPECT_GE(nv, 2); + + // Blend surface should evaluate to finite values + Point3D p = blend.evaluate(0.5, 0.5); + EXPECT_TRUE(std::isfinite(p.x())); + EXPECT_TRUE(std::isfinite(p.y())); + EXPECT_TRUE(std::isfinite(p.z())); +} + +TEST(SurfaceExtensionTest, BlendSurfacesNoCommonBoundary) { + auto plane1 = make_planar_surface(); + std::vector> grid2 = { + {Point3D(3, 0, 0), Point3D(3, 1, 0)}, + {Point3D(4, 0, 0), Point3D(4, 1, 0)} + }; + NurbsSurface plane2(grid2, {0, 0, 1, 1}, {0, 0, 1, 1}, {}, 1, 1); + + auto blend = blend_surfaces(plane1, plane2, 0.3); + auto [nu, nv] = blend.num_control_points(); + // No common boundary → degenerate + EXPECT_TRUE(nu == 0 || nv == 0); +} + +// =========================================================================== +// is_g1_continuous — G1 连续性检查 +// =========================================================================== + +TEST(SurfaceExtensionTest, G1Continuous_AdjacentPlanes) { + auto plane1 = make_planar_surface(); + std::vector> grid2 = { + {Point3D(1, 0, 0), Point3D(1, 1, 0)}, + {Point3D(2, 0, 0), Point3D(2, 1, 0)} + }; + NurbsSurface plane2(grid2, {0, 0, 1, 1}, {0, 0, 1, 1}, {}, 1, 1); + + // They share x=1 boundary: plane1 umax (edge=1) matches plane2 umin (edge=0) + EXPECT_TRUE(is_g1_continuous(plane1, plane2, 1)); +} + +TEST(SurfaceExtensionTest, G1Continuous_NonMatching) { + auto plane1 = make_planar_surface(); + // plane2 is angled → normals won't match + std::vector> grid2 = { + {Point3D(1, 0, 0), Point3D(1, 1, 1)}, + {Point3D(2, 0, 0), Point3D(2, 1, 1)} + }; + NurbsSurface plane2(grid2, {0, 0, 1, 1}, {0, 0, 1, 1}, {}, 1, 1); + + // Normals differ → not G1 + EXPECT_FALSE(is_g1_continuous(plane1, plane2, 1)); +}