From 66777e08399f2c9a7da4ba9ee3eafac4334a2ac7 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 22:08:38 +0800 Subject: [PATCH] feat(v6.1): 5-axis CAM + reverse engineering + FEA mesh generation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit v6.1.1 — 5-Axis CAM: - cam_5axis.h/.cpp: swarf_machining, multi_axis_roughing/finishing - 4 tool axis strategies, AC/BC/AB machine configs - inverse_kinematics, collision detection (holder+shank) - 30 tests, compilation passes v6.1.2 — Reverse Engineering: - point_cloud.h/.cpp: PointCloud class, PLY/E57/PTX loading - voxel_downsample, statistical_outlier_removal, kNN normal estimation - reverse_engineering.h/.cpp: Poisson surface reconstruction - mesh_to_nurbs_surface (quad remesh + LSQ fitting) - fit_plane (SVD), hole_filling, curvature-aware sampling - 19 tests v6.1.3 — FEA Mesh Generation: - fea_mesh.h/.cpp: tetrahedral (Delaunay 3D), boundary layer (prism) - hexahedral (sweep/extrude), adaptive refinement - MeshQuality: skewness, aspect_ratio, Jacobian, orthogonality - export_abaqus/ansys/nastran formats - 15 tests, 7/8 quality metrics verified --- include/vde/core/cam_5axis.h | 275 +++++++ include/vde/mesh/fea_mesh.h | 409 +++++++++++ include/vde/mesh/point_cloud.h | 127 ++++ include/vde/mesh/reverse_engineering.h | 235 ++++++ src/CMakeLists.txt | 4 + src/core/cam_5axis.cpp | 634 ++++++++++++++++ src/mesh/fea_mesh.cpp | 938 ++++++++++++++++++++++++ src/mesh/point_cloud.cpp | 731 ++++++++++++++++++ src/mesh/reverse_engineering.cpp | 770 +++++++++++++++++++ tests/core/CMakeLists.txt | 1 + tests/core/test_cam_5axis.cpp | 431 +++++++++++ tests/mesh/CMakeLists.txt | 2 + tests/mesh/test_fea_mesh.cpp | 356 +++++++++ tests/mesh/test_reverse_engineering.cpp | 268 +++++++ 14 files changed, 5181 insertions(+) create mode 100644 include/vde/core/cam_5axis.h create mode 100644 include/vde/mesh/fea_mesh.h create mode 100644 include/vde/mesh/point_cloud.h create mode 100644 include/vde/mesh/reverse_engineering.h create mode 100644 src/core/cam_5axis.cpp create mode 100644 src/mesh/fea_mesh.cpp create mode 100644 src/mesh/point_cloud.cpp create mode 100644 src/mesh/reverse_engineering.cpp create mode 100644 tests/core/test_cam_5axis.cpp create mode 100644 tests/mesh/test_fea_mesh.cpp create mode 100644 tests/mesh/test_reverse_engineering.cpp diff --git a/include/vde/core/cam_5axis.h b/include/vde/core/cam_5axis.h new file mode 100644 index 0000000..9118de7 --- /dev/null +++ b/include/vde/core/cam_5axis.h @@ -0,0 +1,275 @@ +#pragma once +/** + * @file cam_5axis.h + * @brief 5 轴联动 CAM 加工策略 + * + * 在 3 轴 cam_strategies.h 的基础上扩展为 5 轴联动： + * - 刀位点 + 刀轴矢量 (CLData) + * - 侧刃加工 (swarf_machining) + * - 多轴粗加工 3+2 定位 (multi_axis_roughing) + * - 多轴精加工 + 碰撞检测 (multi_axis_finishing) + * - 机床运动学逆解 (inverse_kinematics) + * - 刀轴策略枚举 (ToolAxisStrategy) + * + * @ingroup core + */ +#include "vde/core/point.h" +#include "vde/core/aabb.h" +#include "vde/core/cam_toolpath.h" +#include "vde/core/cam_strategies.h" +#include "vde/mesh/halfedge_mesh.h" +#include +#include +#include +#include + +namespace vde::brep { +class BrepModel; +} // namespace vde::brep + +namespace vde::curves { +class NurbsSurface; +} // namespace vde::curves + +namespace vde::core { + +using core::Point3D; +using core::Vector3D; +using core::AABB3D; + +// ═══════════════════════════════════════════════════════════════════════════ +// 5 轴刀位数据 +// ═══════════════════════════════════════════════════════════════════════════ + +/// 单个刀位点 + 刀轴矢量 (Cutter Location Data) +struct CLData { + Point3D position; ///< 刀尖点/TCP 位置 + Vector3D tool_axis; ///< 刀轴单位矢量（从刀尖指向刀柄） + double feed_rate = 500.0; ///< 进给速度 (mm/min) + + /// 便捷构造 + CLData() : position(Point3D::Zero()), tool_axis(Vector3D::UnitZ()) {} + CLData(const Point3D& pos, const Vector3D& axis_, double feed = 500.0) + : position(pos), tool_axis(axis_.normalized()), feed_rate(feed) {} +}; + +/// 5 轴刀路 — CLData 集合 +struct AxisToolpath5 { + std::string name; + std::vector points; ///< 刀位点序列 + double safe_height = 50.0; ///< 安全高度 (mm) + + /// 刀位点数量 + [[nodiscard]] size_t size() const { return points.size(); } + + /// 是否为空 + [[nodiscard]] bool empty() const { return points.empty(); } +}; + +// ═══════════════════════════════════════════════════════════════════════════ +// 刀轴策略 +// ═══════════════════════════════════════════════════════════════════════════ + +/// 刀轴控制策略 +enum class ToolAxisStrategy { + FixedAngle, ///< 固定角度：刀轴保持与 Z 轴固定偏角 + LeadLag, ///< 导前/滞后角：沿切削方向前倾或后倾 + Tilted, ///< 侧倾：垂直于切削方向倾斜 + NormalToSurface ///< 法向：刀轴始终沿曲面法向 +}; + +// ═══════════════════════════════════════════════════════════════════════════ +// 机床配置 +// ═══════════════════════════════════════════════════════════════════════════ + +/// 机床结构类型 +enum class MachineType { + HeadHead, ///< 双摆头：A/B 轴均在主轴侧 + TableTable, ///< 双转台：A/C 轴均在工作台侧 + HeadTable ///< 摆头+转台：主轴侧一轴 + 工作台侧一轴 +}; + +/// 旋转轴配置 +enum class MachineAxisConfig { + AC, ///< A 轴（绕 X 旋转）+ C 轴（绕 Z 旋转） + BC, ///< B 轴（绕 Y 旋转）+ C 轴（绕 Z 旋转） + AB ///< A 轴（绕 X 旋转）+ B 轴（绕 Y 旋转） +}; + +/// 机床运动学配置 +struct MachineConfig { + MachineType machine_type = MachineType::TableTable; + MachineAxisConfig axis_config = MachineAxisConfig::AC; + + /// 旋转轴行程限制（度） + double a_min = -120.0; + double a_max = 120.0; + double b_min = -120.0; + double b_max = 120.0; + double c_min = -9999.0; // C 轴通常无限制（360°连续） + double c_max = 9999.0; + + /// 主轴端到旋转中心的偏移 (mm) + Vector3D pivot_offset = Vector3D::Zero(); + + /// 刀具长度 (mm) + double tool_length = 100.0; + + /// 检查旋转角是否在行程范围内 + [[nodiscard]] bool in_range(double a, double b, double c) const; + + /// 解析刀轴矢量为旋转角（根据 axis_config） + /// @return (a_deg, b_deg, c_deg)，失败返回 std::nullopt + [[nodiscard]] std::optional> + axis_to_angles(const Vector3D& tool_axis) const; +}; + +// ═══════════════════════════════════════════════════════════════════════════ +// 5 轴加工参数 +// ═══════════════════════════════════════════════════════════════════════════ + +/// 侧刃加工参数 +struct SwarfParams { + double step_down = 1.0; ///< 每层切深 (mm) + double step_over = 0.5; ///< 横向步距 (mm) + double feed_rate = 600.0; ///< 进给速度 (mm/min) + double safe_height = 50.0; ///< 安全高度 (mm) + double tilt_angle = 0.0; ///< 刀轴相对于曲面法线的偏转角 (度) +}; + +/// 多轴粗加工参数 +struct MultiAxisRoughingParams { + double step_down = 2.0; ///< 每层切深 (mm) + double step_over = 3.0; ///< 横向步距 (mm) + double stock_to_leave = 1.0; ///< 留量 (mm) + double feed_rate = 800.0; ///< 进给速度 (mm/min) + double safe_height = 50.0; ///< 安全高度 (mm) + /// 3+2 定位角度 — 固定刀轴方向（相对于 Z 轴的倾斜角, 度） + double tilt_angle = 15.0; + double lead_angle = 0.0; ///< 导前角 (度) +}; + +/// 多轴精加工参数 +struct MultiAxisFinishingParams { + double step_over = 0.3; ///< 行距 (mm) + double feed_rate = 500.0; ///< 进给速度 (mm/min) + double safe_height = 50.0; ///< 安全高度 (mm) + ToolAxisStrategy axis_strategy = ToolAxisStrategy::NormalToSurface; + double max_tilt = 30.0; ///< 最大倾斜角 (度) + double lead_angle = 3.0; ///< 导前角 (度) + /// 碰撞检测参数 + bool collision_check = true; ///< 是否启用碰撞检测 + double shank_clearance = 2.0; ///< 刀柄安全间隙 (mm) + double holder_clearance = 5.0; ///< 刀杆安全间隙 (mm) +}; + +// ═══════════════════════════════════════════════════════════════════════════ +// 5 轴加工策略函数 +// ═══════════════════════════════════════════════════════════════════════════ + +/// 侧刃加工 (Swarf Machining) +/// +/// 使用刀具侧刃贴合曲面进行切削，刀轴沿曲面法线偏转。 +/// 适用于直纹面侧壁加工。 +/// +/// @param surface 目标 NURBS 曲面 +/// @param tool 刀具参数 +/// @param params 侧刃加工参数 +/// @return 5 轴刀路 +/// +/// @ingroup core +[[nodiscard]] AxisToolpath5 swarf_machining( + const curves::NurbsSurface& surface, + const Tool& tool, + const SwarfParams& params); + +/// 多轴粗加工 (3+2 定位加工) +/// +/// 使用 3+2 定位方式对实体进行分层粗加工。 +/// 自动优化刀轴方向以最大化材料去除率，同时避免干涉。 +/// +/// @param body B-Rep 实体 +/// @param tool 刀具参数 +/// @param params 多轴粗加工参数 +/// @param config 机床运动学配置（用于刀轴方向可行性检查） +/// @return 5 轴刀路 +/// +/// @ingroup core +[[nodiscard]] AxisToolpath5 multi_axis_roughing( + const brep::BrepModel& body, + const Tool& tool, + const MultiAxisRoughingParams& params, + const MachineConfig& config = MachineConfig{}); + +/// 多轴精加工 +/// +/// 对曲面进行 5 轴联动精加工。刀轴平滑变化， +/// 包含刀柄和刀杆的碰撞检测。 +/// +/// @param surface 目标 NURBS 曲面 +/// @param body 完整 B-Rep 模型（碰撞检测用） +/// @param tool 刀具参数 +/// @param params 多轴精加工参数 +/// @param config 机床运动学配置 +/// @return 5 轴刀路 +/// +/// @ingroup core +[[nodiscard]] AxisToolpath5 multi_axis_finishing( + const curves::NurbsSurface& surface, + const brep::BrepModel& body, + const Tool& tool, + const MultiAxisFinishingParams& params, + const MachineConfig& config = MachineConfig{}); + +/// 机床运动学逆解 +/// +/// 将工件坐标系下的刀路 (CLData) 转换为机床各轴坐标。 +/// 根据机床配置 (AC/BC/AB + 双转台/摆头) 计算 X/Y/Z/A/B/C。 +/// +/// @param toolpath 工件坐标系下的 5 轴刀路 +/// @param config 机床运动学配置 +/// @return 机床坐标刀路（points 中 position 存 X/Y/Z，tool_axis 存 A/B/C 度） +/// +/// @ingroup core +[[nodiscard]] AxisToolpath5 inverse_kinematics( + const AxisToolpath5& toolpath, + const MachineConfig& config); + +/// 刀轴碰撞检测 +/// +/// 检查刀柄和刀杆是否与工件发生碰撞。 +/// +/// @param position 刀尖位置 +/// @param tool_axis 刀轴方向 +/// @param body B-Rep 工件 +/// @param tool 刀具参数 +/// @param shank_clearance 刀柄安全间隙 +/// @param holder_clearance 刀杆安全间隙 +/// @return true 如果发生碰撞 +/// +/// @ingroup core +[[nodiscard]] bool check_collision( + const Point3D& position, + const Vector3D& tool_axis, + const brep::BrepModel& body, + const Tool& tool, + double shank_clearance = 2.0, + double holder_clearance = 5.0); + +/// 刀轴方向优化 +/// +/// 给定候选刀轴方向列表，选择最优方向： +/// - 优先选择更接近 Z 轴（减小旋转轴行程） +/// - 避免超出机床行程的配置 +/// - 考虑切削条件（避免刀尖切削） +/// +/// @param candidates 候选刀轴方向 +/// @param config 机床配置 +/// @return 最优刀轴方向（归一化），若无可行方向返回 Z 轴 +/// +/// @ingroup core +[[nodiscard]] Vector3D optimize_tool_axis( + const std::vector& candidates, + const MachineConfig& config); + +} // namespace vde::core diff --git a/include/vde/mesh/fea_mesh.h b/include/vde/mesh/fea_mesh.h new file mode 100644 index 0000000..727701a --- /dev/null +++ b/include/vde/mesh/fea_mesh.h @@ -0,0 +1,409 @@ +#pragma once +/** + * @file fea_mesh.h + * @brief 高质量FEA有限元网格生成模块 + * + * 提供四面体、边界层棱柱、六面体网格生成、自适应细化、 + * 网格质量评估及CAE格式（Abaqus/Ansys/Nastran）导出。 + * + * ## 功能概览 + * - tetrahedral_mesh: Delaunay 3D + 质量优化四面体网格 + * - boundary_layer_mesh: 沿壁面法向拉伸棱柱边界层 + * - hexahedral_mesh: sweep/extrude 六面体网格 + * - adaptive_refinement: 基于误差估计的局部细化 + * - fea_quality_report: skewness/aspect_ratio/jacobian/orthogonality 报告 + * - export_abaqus / export_ansys / export_nastran: CAE 格式导出 + * + * @ingroup mesh + */ +#include "vde/mesh/halfedge_mesh.h" +#include "vde/mesh/delaunay_3d.h" +#include "vde/brep/brep.h" +#include +#include +#include +#include +#include + +namespace vde::mesh { +using core::Point3D; +using core::Vector3D; + +// ════════════════════════════════════════════════════════ +// FEA 单元类型 +// ════════════════════════════════════════════════════════ + +/** + * @brief FEA 单元拓扑类型 + * @ingroup mesh + */ +enum class FEAElementType { + Tet4, ///< 4节点线性四面体 (C3D4) + Tet10, ///< 10节点二次四面体 (C3D10) + Hex8, ///< 8节点线性六面体 (C3D8) + Hex20, ///< 20节点二次六面体 (C3D20) + Wedge6, ///< 6节点线性三棱柱/楔形 (C3D6) + Wedge15, ///< 15节点二次三棱柱 (C3D15) +}; + +/// 获取单元类型对应的每单元节点数 +inline int nodes_per_element(FEAElementType t) { + switch (t) { + case FEAElementType::Tet4: return 4; + case FEAElementType::Tet10: return 10; + case FEAElementType::Hex8: return 8; + case FEAElementType::Hex20: return 20; + case FEAElementType::Wedge6: return 6; + case FEAElementType::Wedge15:return 15; + } + return 0; +} + +/// 获取单元类型名称 +inline const char* element_type_name(FEAElementType t) { + switch (t) { + case FEAElementType::Tet4: return "Tet4"; + case FEAElementType::Tet10: return "Tet10"; + case FEAElementType::Hex8: return "Hex8"; + case FEAElementType::Hex20: return "Hex20"; + case FEAElementType::Wedge6: return "Wedge6"; + case FEAElementType::Wedge15:return "Wedge15"; + } + return "Unknown"; +} + +// ════════════════════════════════════════════════════════ +// FEAMesh 数据结构 +// ════════════════════════════════════════════════════════ + +/** + * @brief FEA有限元网格 + * + * 存储顶点、单元连通性、边界面及单元类型信息。 + * 支持Tet4/Tet10/Hex8/Wedge6等常见FEA单元。 + * + * @ingroup mesh + */ +struct FEAMesh { + std::vector vertices; ///< 顶点坐标数组 + std::vector> elements; ///< 各单元顶点索引（长度 = nodes_per_element） + FEAElementType element_type = FEAElementType::Tet4; ///< 单元类型 + std::vector material_ids; ///< 各单元材料ID（可选，默认0） + std::vector> boundary_faces;///< 边界面（用于施加BC），顶点索引数组 + std::vector boundary_face_tags; ///< 边界面标签（与boundary_faces一一对应） + + /** @brief 单元总数 */ + [[nodiscard]] size_t num_elements() const { return elements.size(); } + /** @brief 顶点总数 */ + [[nodiscard]] size_t num_vertices() const { return vertices.size(); } + /** @brief 边界面总数 */ + [[nodiscard]] size_t num_boundary_faces() const { return boundary_faces.size(); } + /** @brief 每单元节点数 */ + [[nodiscard]] int npe() const { return nodes_per_element(element_type); } +}; + +// ════════════════════════════════════════════════════════ +// 网格生成参数 +// ════════════════════════════════════════════════════════ + +/** + * @brief 四面体网格生成参数 + * @ingroup mesh + */ +struct TetMeshParams { + double min_angle = 15.0; ///< 最小二面角（度），太小导致刚度矩阵病态 + double max_aspect_ratio = 5.0; ///< 最大长宽比（外接球半径/最短边） + double max_size = 0.0; ///< 最大单元尺寸，0 = 自动（包围盒对角线的 5%） + double min_size = 0.0; ///< 最小单元尺寸，0 = max_size/10 + int quality_iterations = 3; ///< 质量优化迭代次数（Laplacian光顺+边翻转） + bool keep_surface = true; ///< 是否保持原表面顶点位置 +}; + +/** + * @brief 边界层网格参数 + * @ingroup mesh + */ +struct BLPParams { + double first_cell_height = 0.001; ///< 第一层棱柱高度（壁面法向距离） + double growth_rate = 1.2; ///< 相邻层高度增长率（>1.0） + int num_layers = 5; ///< 棱柱层数 + double total_thickness = 0.0; ///< 边界层总厚度，0 = 自动（first * (growth^n - 1)/(growth - 1)） + bool smooth_transition = true; ///< 是否平滑过渡到内部网格 +}; + +/** + * @brief 六面体网格参数 + * @ingroup mesh + */ +struct HexMeshParams { + int sweep_layers = 1; ///< sweep 方向层数 + double sweep_distance = 0.0; ///< sweep 总距离，0 = 自动（包围盒最小维度） + bool use_quad_dominant = true; ///< 源面是否使用四边形主导网格 + double cell_size = 0.0; ///< 目标单元尺寸，0 = 自动 +}; + +// ════════════════════════════════════════════════════════ +// 网格质量评估 +// ════════════════════════════════════════════════════════ + +/** + * @brief FEA网格质量报告 + * + * 包含 skewness（偏斜度）、aspect_ratio（长宽比）、 + * jacobian（雅可比）、orthogonality（正交性）四项 + * 核心FEA网格质量指标。 + * + * @ingroup mesh + */ +struct FEAQualityReport { + /// 偏斜度 (skewness): 实际单元与理想单元的偏差，[0,1]，0=最优 + double min_skewness = 0.0; + double max_skewness = 0.0; + double avg_skewness = 0.0; + + /// 长宽比 (aspect ratio): 最长边/最短边，(Fluent: 1=最佳) + double min_aspect_ratio = 0.0; + double max_aspect_ratio = 0.0; + double avg_aspect_ratio = 0.0; + + /// 归一化雅可比 (scaled jacobian): [0,1]，1=正单元，0=退化，<0=反转 + double min_jacobian = 0.0; + double max_jacobian = 0.0; + double avg_jacobian = 0.0; + + /// 正交性 (orthogonality): 面法向与质心连线的夹角余弦，[0,1]，1=完美正交 + double min_orthogonality = 0.0; + double max_orthogonality = 0.0; + double avg_orthogonality = 0.0; + + size_t total_elements = 0; ///< 单元总数 + size_t degenerate_elements = 0; ///< 退化单元数 (Jacobian < 0.01) + size_t inverted_elements = 0; ///< 反转单元数 (Jacobian < 0) + + /// 分级统计 (A=优秀 B=良好 C=合格 D=差 F=不合格，按Jacobian) + size_t grade_A = 0, grade_B = 0, grade_C = 0, grade_D = 0, grade_F = 0; +}; + +// ════════════════════════════════════════════════════════ +// 自适应细化 +// ════════════════════════════════════════════════════════ + +/** + * @brief 误差估计函数类型 + * + * 对给定单元返回估计误差值（≥0），值越大表示该单元越需要细化。 + * + * @param element_idx 单元索引 + * @param mesh 所属网格 + * @return 误差估计值（非负） + */ +using ErrorEstimator = std::function; + +/** + * @brief 自适应细化参数 + * @ingroup mesh + */ +struct RefinementParams { + double error_threshold = 0.5; ///< 相对误差阈值（高于此值的单元被细化） + double refinement_fraction = 0.2; ///< 每次细化单元比例（0~1） + int max_iterations = 3; ///< 最大细化迭代次数 + bool smooth_after_refine = true; ///< 细化后是否光顺 +}; + +// ════════════════════════════════════════════════════════ +// 网格生成 API +// ════════════════════════════════════════════════════════ + +/** + * @brief 生成四面体FEA网格（Delaunay 3D + 质量优化） + * + * 对 BrepModel 进行 tessellation 获取表面三角网格， + * 采样内部点后进行 Delaunay 3D 四面体化， + * 再进行 Laplacian 光顺 + 边翻转优化质量。 + * + * @param body B-Rep实体模型 + * @param params 四面体化参数（最小角、长宽比约束等） + * @return FEAMesh 四面体网格 + * + * @code{.cpp} + * auto box = vde::brep::make_box(1, 1, 1); + * TetMeshParams p; + * p.min_angle = 18.0; + * p.max_aspect_ratio = 4.0; + * auto mesh = tetrahedral_mesh(box, p); + * @endcode + * @ingroup mesh + */ +FEAMesh tetrahedral_mesh(const brep::BrepModel& body, + const TetMeshParams& params = {}); + +/** + * @brief 生成边界层棱柱网格 + * + * 沿B-Rep模型壁面法向拉伸多层棱柱单元， + * 用于CFD边界层捕捉（y+控制）。 + * 先对表面进行三角化，然后每层沿法向偏移。 + * + * @param body B-Rep实体模型 + * @param params 边界层参数（首层高度、增长率、层数） + * @return FEAMesh 棱柱网格 (Wedge6) + * + * @code{.cpp} + * BLPParams p; + * p.first_cell_height = 0.0001; + * p.growth_rate = 1.15; + * p.num_layers = 10; + * auto bl = boundary_layer_mesh(wing, p); + * @endcode + * @ingroup mesh + */ +FEAMesh boundary_layer_mesh(const brep::BrepModel& body, + const BLPParams& params = {}); + +/** + * @brief 生成六面体网格（sweep/extrude策略） + * + * 将B-Rep模型的一个面（源面）四边形化， + * 然后沿法向 sweep 生成六面体层。 + * 适用于拉伸体/棱柱体等 sweepable 几何。 + * + * @param body B-Rep实体模型 + * @param params 六面体化参数 + * @return FEAMesh 六面体网格 (Hex8) + * + * @code{.cpp} + * auto cylinder = vde::brep::make_cylinder(1, 3); + * HexMeshParams p; + * p.sweep_layers = 5; + * auto mesh = hexahedral_mesh(cylinder, p); + * @endcode + * @ingroup mesh + */ +FEAMesh hexahedral_mesh(const brep::BrepModel& body, + const HexMeshParams& params = {}); + +/** + * @brief 自适应网格细化 + * + * 基于误差估计函数对高误差区域进行局部细分。 + * 支持边中点分裂细化（适用于Tet4→4×Tet4子单元）。 + * + * @param mesh 输入FEA网格 + * @param estimator 误差估计函数 + * @param params 细化参数 + * @return FEAMesh 细化后的网格 + * + * @code{.cpp} + * auto estimator = [](int ei, const FEAMesh& m) -> double { + * // 基于应力梯度的误差估计 + * return stress_gradient_norm(ei, m); + * }; + * auto refined = adaptive_refinement(mesh, estimator); + * @endcode + * @ingroup mesh + */ +FEAMesh adaptive_refinement(const FEAMesh& mesh, + const ErrorEstimator& estimator, + const RefinementParams& params = {}); + +// ════════════════════════════════════════════════════════ +// 质量评估 API +// ════════════════════════════════════════════════════════ + +/** + * @brief 评估FEA网格质量 + * + * 计算四项核心质量指标： + * - skewness: 偏斜度（实际vs理想单元的偏差） + * - aspect_ratio: 长宽比 + * - jacobian: 归一化雅可比行列式 + * - orthogonality: 正交性（面法向与质心连线夹角） + * + * 自动根据 element_type 选择对应算法（Tet/Hex/Wedge）。 + * + * @param mesh FEA网格 + * @return FEAQualityReport 质量报告 + * + * @ingroup mesh + */ +FEAQualityReport fea_quality_report(const FEAMesh& mesh); + +// ── 单项质量指标 ── + +/** + * @brief 计算单元的偏斜度 + * @param verts 单元顶点数组 + * @param type 单元类型 + * @return skewness [0,1]，0=理想单元 + */ +double element_skewness(const std::vector& verts, FEAElementType type); + +/** + * @brief 计算单元的长宽比 + * @param verts 单元顶点数组 + * @param type 单元类型 + * @return aspect_ratio ≥ 1，1=理想 + */ +double element_aspect_ratio(const std::vector& verts, FEAElementType type); + +/** + * @brief 计算单元的归一化雅可比 + * @param verts 单元顶点数组 + * @param type 单元类型 + * @return scaled_jacobian [0,1]，1=正四面体/正方体，0=退化 + */ +double element_scaled_jacobian(const std::vector& verts, FEAElementType type); + +/** + * @brief 计算单元的正交性 + * @param verts 单元顶点数组 + * @param type 单元类型 + * @return orthogonality [0,1]，1=完全正交（面法向与对面质心连线平行） + */ +double element_orthogonality(const std::vector& verts, FEAElementType type); + +// ════════════════════════════════════════════════════════ +// CAE 格式导出 +// ════════════════════════════════════════════════════════ + +/** + * @brief 导出 Abaqus 输入文件 (.inp) + * + * 生成包含 *NODE 和 *ELEMENT 卡的 Abaqus 输入文件。 + * 支持 Tet4 (C3D4), Tet10 (C3D10), Hex8 (C3D8), Wedge6 (C3D6)。 + * + * @param mesh FEA网格 + * @param filepath 输出文件路径 + * @return 成功返回 true + * + * @ingroup mesh + */ +bool export_abaqus(const FEAMesh& mesh, const std::string& filepath); + +/** + * @brief 导出 Ansys CDB 文件 (.cdb) + * + * 生成包含 NBLOCK 和 EBLOCK 的 Ansys CDB 格式文件。 + * + * @param mesh FEA网格 + * @param filepath 输出文件路径 + * @return 成功返回 true + * + * @ingroup mesh + */ +bool export_ansys(const FEAMesh& mesh, const std::string& filepath); + +/** + * @brief 导出 Nastran 批量数据文件 (.bdf/.dat) + * + * 生成包含 GRID 和 CTETRA/CHEXA/CPENTA 卡的 Nastran 输入文件。 + * + * @param mesh FEA网格 + * @param filepath 输出文件路径 + * @return 成功返回 true + * + * @ingroup mesh + */ +bool export_nastran(const FEAMesh& mesh, const std::string& filepath); + +} // namespace vde::mesh diff --git a/include/vde/mesh/point_cloud.h b/include/vde/mesh/point_cloud.h new file mode 100644 index 0000000..1dddc37 --- /dev/null +++ b/include/vde/mesh/point_cloud.h @@ -0,0 +1,127 @@ +#pragma once +#include "vde/core/point.h" +#include "vde/core/aabb.h" +#include +#include +#include +#include + +namespace vde::mesh { +using core::Point3D; +using core::Vector3D; +using core::AABB3D; + +// ─────────────────────────────────────────────────────────── +// PointCloud 点云类 +// ─────────────────────────────────────────────────────────── + +/** + * @brief 点云数据结构 + * + * 存储顶点位置、法线、颜色（可选）及包围盒信息。 + * 支持 PLY/E57/PTX 加载，体素降采样、统计离群点剔除、kNN 法线估计。 + * + * @ingroup mesh + */ +class PointCloud { +public: + std::vector vertices; ///< 顶点位置 + std::vector normals; ///< 顶点法线（归一化，空数组表示未计算） + std::vector> colors; ///< RGB 颜色（0-255，空数组表示无颜色） + + PointCloud() = default; + + /// 从顶点列表构造（无颜色/法线） + explicit PointCloud(std::vector pts) : vertices(std::move(pts)) {} + + /// 顶点数 + [[nodiscard]] size_t size() const { return vertices.size(); } + /// 是否有法线 + [[nodiscard]] bool has_normals() const { return !normals.empty(); } + /// 是否有颜色 + [[nodiscard]] bool has_colors() const { return !colors.empty(); } + + /// 包围盒 + [[nodiscard]] AABB3D bounds() const; + + // ── I/O ────────────────────────────────────────────── + + /// 加载 PLY 点云文件（ASCII 或 binary little-endian） + static PointCloud load_ply(const std::string& path); + + /// 加载 E57 点云文件（简化解析：仅提取 CartesianBounds + points3D） + static PointCloud load_e57(const std::string& path); + + /// 加载 PTX 格式点云（Leica 扫描仪原始格式） + static PointCloud load_ptx(const std::string& path); + + /// 保存为 PLY 文件 + bool save_ply(const std::string& path) const; + + // ── 滤波 / 降采样 ──────────────────────────────────── + + /** + * @brief 体素降采样 + * @param voxel_size 体素边长 + * @return 降采样后的点云（每个体素用质心代替） + */ + [[nodiscard]] PointCloud voxel_downsample(double voxel_size) const; + + /** + * @brief 统计离群点剔除（SOR） + * @param k 邻域点数（默认 6） + * @param std_ratio 标准差倍率阈值（默认 1.0），超过 mean+k*std_ratio 的点标记为离群 + * @return 剔除离群点后的点云 + */ + [[nodiscard]] PointCloud statistical_outlier_removal(int k = 6, double std_ratio = 1.0) const; + + // ── 法线估计 ────────────────────────────────────────── + + /** + * @brief kNN 法线估计 + * + * 对每个点，找到 k 个最近邻点，用 PCA 估计法线（最小特征值对应方向）。 + * 法线方向通过视角向量（点→原点方向）统一朝向。 + * + * @param k 最近邻点数（≥ 3） + * @return 估计了法线的新点云 + */ + [[nodiscard]] PointCloud normal_estimation(int k = 10) const; + + // ── 曲率感知采样 ────────────────────────────────────── + + /** + * @brief 曲率感知降采样 + * + * 计算每点局部曲率近似，按曲率加权采样：高曲率区域保留更多点。 + * + * @param target_count 目标点数 + * @return 采样后的索引数组 + */ + [[nodiscard]] std::vector curvature_aware_sampling(size_t target_count) const; + + // ── 平滑滤波 ────────────────────────────────────────── + + /** + * @brief 双边滤波平滑 + * + * 对点云位置进行法向加权的双边滤波去噪，保持尖锐特征。 + * + * @param sigma_c 空间高斯 σ（0 = 自动按平均点距计算） + * @param sigma_n 法向高斯 σ（弧度） + * @param iterations 迭代次数 + * @return 平滑后的新点云 + */ + [[nodiscard]] PointCloud smoothing_filter(double sigma_c = 0.0, + double sigma_n = 0.3, + int iterations = 3) const; + +private: + /** 构建 k-d 树索引，返回每个顶点的 k 最近邻索引 */ + [[nodiscard]] std::vector> k_nearest_neighbors(size_t k) const; + + /** 协方差矩阵 PCA：计算三个特征值，返回最小特征值对应的特征向量 */ + [[nodiscard]] static Vector3D pca_normal(const std::vector& local_pts); +}; + +} // namespace vde::mesh diff --git a/include/vde/mesh/reverse_engineering.h b/include/vde/mesh/reverse_engineering.h new file mode 100644 index 0000000..4aee24b --- /dev/null +++ b/include/vde/mesh/reverse_engineering.h @@ -0,0 +1,235 @@ +#pragma once +#include "vde/mesh/halfedge_mesh.h" +#include "vde/mesh/point_cloud.h" +#include "vde/curves/nurbs_surface.h" +#include +#include + +namespace vde::mesh { +using core::Point3D; +using core::Vector3D; +using core::AABB3D; +using curves::NurbsSurface; + +// ─────────────────────────────────────────────────────────── +// 参数结构体 +// ─────────────────────────────────────────────────────────── + +/** + * @brief Poisson 曲面重建参数 + * @ingroup mesh + */ +struct PoissonReconParams { + int max_depth = 8; ///< 八叉树最大深度（越高越精细） + int min_depth = 3; ///< 八叉树最小深度 + double point_weight = 4.0; ///< 插值点权重（越高越贴近输入点） + int solver_divide = 8; ///< 求解器分块大小 + int iso_divide = 8; ///< 等值面提取分块大小 + double samples_per_node = 1.5; ///< 每节点采样数 + bool scale_fit = true; ///< 是否缩放到单位立方体 + double confidence = 0.0; ///< 点置信度（0=均匀权重） +}; + +/** + * @brief 四边形重网格化参数 + * @ingroup mesh + */ +struct QuadRemeshParams { + int target_quads = 400; ///< 目标四边形数量 + double edge_length = 0.0; ///< 目标边长（0=自动） + double regularity = 0.5; ///< 规则性权重 [0,1] + double sharpness = 0.5; ///< 特征保持权重 [0,1] + int iterations = 10; ///< 迭代次数 + bool adaptive_size = false; ///< 自适应边长（曲率区域更密） +}; + +/** + * @brief NURBS 拟合参数 + * @ingroup mesh + */ +struct NURBSFitParams { + int degree_u = 3; ///< u 向阶次 + int degree_v = 3; ///< v 向阶次 + int num_cp_u = 8; ///< u 向控制点数（≥ degree_u+1） + int num_cp_v = 8; ///< v 向控制点数（≥ degree_v+1） + double smoothness = 0.001; ///< 平滑项权重（Tikhonov 正则化） + int max_iterations = 50; ///< 最小二乘迭代上限 + double tolerance = 1e-6; ///< 收敛容差 + bool use_weights = true; ///< 是否使用有理权重拟合（NURBS vs B-spline） +}; + +/** + * @brief 点云→网格综合参数 + * @ingroup mesh + */ +struct PointCloudToMeshParams { + int knn = 10; ///< 法线估计 kNN + double voxel_size = 0.0; ///< 预降采样体素尺寸（0=跳过） + bool remove_outliers = true; ///< 是否先剔除统计离群点 + PoissonReconParams poisson; ///< Poisson 重建参数 +}; + +/** + * @brief 网格→NURBS 综合参数 + * @ingroup mesh + */ +struct MeshToNURBSParams { + QuadRemeshParams quad; ///< 四边形重网格化参数 + NURBSFitParams nurbs; ///< NURBS 拟合参数 + bool smooth_before_remesh = false; ///< 重网格化前是否平滑 +}; + +// ─────────────────────────────────────────────────────────── +// 拟合结果 +// ─────────────────────────────────────────────────────────── + +/// 平面拟合结果 +struct PlaneFitResult { + Point3D origin; ///< 平面上一点（点云质心） + Vector3D normal; ///< 平面法向量（单位向量） + double rms_error; ///< RMS 拟合误差 + bool valid; ///< 拟合是否有效 +}; + +/// 曲面片拟合结果 +struct PatchFitResult { + NurbsSurface surface; ///< 拟合的 NURBS 曲面 + double rms_error; ///< RMS 拟合误差 + double max_error; ///< 最大偏差 + bool valid; ///< 拟合是否有效 +}; + +// ─────────────────────────────────────────────────────────── +// 核心函数声明 +// ─────────────────────────────────────────────────────────── + +/** + * @brief 点云 → 三角网格（Poisson 曲面重建） + * + * 完整流水线：法线估计 → Poisson 隐式曲面重建 → 等值面提取。 + * + * @param cloud 输入点云（可无法线，内部自动估计） + * @param params 控制参数 + * @return 重建的三角网格 + * + * @code{.cpp} + * PointCloudToMeshParams p; + * p.knn = 10; p.voxel_size = 0.02; + * auto mesh = point_cloud_to_mesh(cloud, p); + * @endcode + * @ingroup mesh + */ +HalfedgeMesh point_cloud_to_mesh(const PointCloud& cloud, + const PointCloudToMeshParams& params = {}); + +/** + * @brief 网格 → NURBS 曲面 + * + * 流水线：四边形重网格化 → 参数化 → 控制点最小二乘拟合。 + * + * @param mesh 输入三角网格 + * @param params 控制参数 + * @return 拟合的 NURBS 曲面 + * + * @code{.cpp} + * MeshToNURBSParams p; + * p.quad.target_quads = 200; + * p.nurbs.num_cp_u = 6; p.nurbs.num_cp_v = 6; + * auto surf = mesh_to_nurbs_surface(mesh, p); + * @endcode + * @ingroup mesh + */ +NurbsSurface mesh_to_nurbs_surface(const HalfedgeMesh& mesh, + const MeshToNURBSParams& params = {}); + +/** + * @brief 点云平面拟合（SVD 最小二乘） + * + * @param points 输入点列表 + * @return 平面拟合结果（origin, normal, rms_error） + * + * @code{.cpp} + * auto plane = fit_plane(points); + * // plane.normal 是单位法向量，plane.origin 是质心 + * @endcode + * @ingroup mesh + */ +PlaneFitResult fit_plane(const std::vector& points); + +/** + * @brief 点云曲面片 NURBS 拟合 + * + * 将点云拟合为指定阶次的 NURBS 曲面片：先将点投影到拟合平面上， + * 在参数域中均匀采样构造控制网格，再用最小二乘优化控制点位置。 + * + * @param points 输入点列表 + * @param degree_u u 向阶次 + * @param degree_v v 向阶次 + * @param num_cp_u u 向控制点数 + * @param num_cp_v v 向控制点数 + * @return 拟合结果（NURBS 曲面 + 误差） + * + * @ingroup mesh + */ +PatchFitResult fit_patch(const std::vector& points, + int degree_u = 3, int degree_v = 3, + int num_cp_u = 8, int num_cp_v = 8); + +/** + * @brief 点云曲率感知采样 + * + * 计算局部曲率近似，高曲率区域保留更多点，低曲率区域用较少点表示。 + * + * @param points 输入点列表 + * @param target_count 目标保留点数 + * @return 采样后的点索引数组 + * + * @ingroup mesh + */ +std::vector curvature_aware_sampling(const std::vector& points, + size_t target_count); + +/** + * @brief 网格孔洞填充 + * + * 检测边界环，对每个孔洞用 fan 三角化（或最小面积三角化）填充。 + * + * @param mesh 输入网格 + * @param max_hole_edges 最大孔洞边界边数（超过此数的孔洞不填充，0=无限制） + * @return 填充后的网格 + * + * @ingroup mesh + */ +HalfedgeMesh hole_filling(const HalfedgeMesh& mesh, int max_hole_edges = 0); + +/** + * @brief 点云平滑滤波 + * + * 对点云位置进行各向同性拉普拉斯 + 双边滤波去噪。 + * + * @param points 输入点列表 + * @param sigma 高斯核 σ（控制平滑强度） + * @param iterations 迭代次数 + * @return 平滑后的点列表 + * + * @ingroup mesh + */ +std::vector smoothing_filter(const std::vector& points, + double sigma = 0.5, int iterations = 3); + +// ─────────────────────────────────────────────────────────── +// 法线估计（独立函数） +// ─────────────────────────────────────────────────────────── + +/** + * @brief 点云法线估计（kNN + PCA） + * + * @param points 输入点列表 + * @param k 最近邻点数 + * @return 估计的法线向量列表 + * + * @ingroup mesh + */ +std::vector estimate_normals(const std::vector& points, int k = 10); + +} // namespace vde::mesh diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 59e00e1..9f62861 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -29,6 +29,7 @@ add_library(vde_core STATIC core/icp.cpp core/cam_toolpath.cpp core/cam_strategies.cpp + core/cam_5axis.cpp core/exact_predicates.cpp ) target_include_directories(vde_core @@ -83,6 +84,9 @@ add_library(vde_mesh STATIC mesh/marching_cubes.cpp mesh/mesh_lod.cpp mesh/parallel_mc.cpp + mesh/reverse_engineering.cpp + mesh/point_cloud.cpp + mesh/fea_mesh.cpp ) target_include_directories(vde_mesh PUBLIC ${CMAKE_SOURCE_DIR}/include diff --git a/src/core/cam_5axis.cpp b/src/core/cam_5axis.cpp new file mode 100644 index 0000000..8a1120c --- /dev/null +++ b/src/core/cam_5axis.cpp @@ -0,0 +1,634 @@ +#include "vde/core/cam_5axis.h" +#include "vde/brep/brep.h" +#include "vde/curves/nurbs_surface.h" +#include "vde/brep/modeling.h" +#include "vde/core/aabb.h" +#include +#include +#include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +namespace vde::core { + +// ═══════════════════════════════════════════════════════════════════════════ +// 内部辅助函数 +// ═══════════════════════════════════════════════════════════════════════════ + +namespace { + +inline double deg2rad(double deg) { return deg * M_PI / 180.0; } +inline double rad2deg(double rad) { return rad * 180.0 / M_PI; } + +/// 将角度归一化到 [-180, 180] 范围 +double normalize_angle(double deg) { + while (deg > 180.0) deg -= 360.0; + while (deg < -180.0) deg += 360.0; + return deg; +} + +/// 限制值在 [lo, hi] 区间 +double clamp(double v, double lo, double hi) { + return std::max(lo, std::min(hi, v)); +} + +/// 计算 NURBS 曲面的参数域边界 +void surface_domain(const curves::NurbsSurface& surf, + double& u0, double& u1, double& v0, double& v1) { + const auto& ku = surf.knots_u(); + const auto& kv = surf.knots_v(); + int pu = surf.degree_u(); + int pv = surf.degree_v(); + u0 = ku[pu]; + u1 = ku[ku.size() - 1 - pu]; + v0 = kv[pv]; + v1 = kv[kv.size() - 1 - pv]; +} + +/// 在曲面上采样点列表 +std::vector sample_surface_grid( + const curves::NurbsSurface& surf, int nu, int nv) { + std::vector pts; + pts.reserve(nu * nv); + double u0, u1, v0, v1; + surface_domain(surf, u0, u1, v0, v1); + for (int i = 0; i < nu; ++i) { + double u = u0 + (u1 - u0) * static_cast(i) / static_cast(nu - 1); + for (int j = 0; j < nv; ++j) { + double v = v0 + (v1 - v0) * static_cast(j) / static_cast(nv - 1); + pts.push_back(surf.evaluate(u, v)); + } + } + return pts; +} + +/// 射线-三角形求交 (Möller-Trumbore) +bool ray_triangle_intersect(const Point3D& origin, + const Vector3D& dir, + const Point3D& v0, + const Point3D& v1, + const Point3D& v2, + double& t_out) { + Vector3D e1 = v1 - v0; + Vector3D e2 = v2 - v0; + Vector3D pvec = dir.cross(e2); + double det = e1.dot(pvec); + if (std::abs(det) < 1e-12) return false; + double inv_det = 1.0 / det; + Vector3D tvec = origin - v0; + double u = tvec.dot(pvec) * inv_det; + if (u < 0.0 || u > 1.0) return false; + Vector3D qvec = tvec.cross(e1); + double v = dir.dot(qvec) * inv_det; + if (v < 0.0 || u + v > 1.0) return false; + t_out = e2.dot(qvec) * inv_det; + return t_out >= 0.0; +} + +/// 射线-网格求交（最近交点） +bool ray_mesh_intersect(const Point3D& origin, + const Vector3D& dir, + const mesh::HalfedgeMesh& mesh, + double max_dist, + double& t_out) { + t_out = max_dist; + bool hit = false; + for (size_t fi = 0; fi < mesh.num_faces(); ++fi) { + auto fv = mesh.face_vertices(static_cast(fi)); + if (fv.size() < 3) continue; + Point3D p0 = mesh.vertex(fv[0]); + Point3D p1 = mesh.vertex(fv[1]); + Point3D p2 = mesh.vertex(fv[2]); + double t; + if (ray_triangle_intersect(origin, dir, p0, p1, p2, t) && t < t_out) { + t_out = t; + hit = true; + } + } + return hit; +} + +/// 计算工件 Z 范围 +std::pair mesh_z_range(const mesh::HalfedgeMesh& mesh) { + double z_min = std::numeric_limits::max(); + double z_max = std::numeric_limits::lowest(); + for (size_t vi = 0; vi < mesh.num_vertices(); ++vi) { + auto v = mesh.vertex(static_cast(vi)); + z_min = std::min(z_min, v.z()); + z_max = std::max(z_max, v.z()); + } + return {z_min, z_max}; +} + +/// 生成切触点序列（沿参数线扫描，返回 (position, normal) 对） +std::vector> generate_contact_points( + const curves::NurbsSurface& surface, + double step_over, + int samples_per_pass = 20) { + std::vector> result; + double u0, u1, v0, v1; + surface_domain(surface, u0, u1, v0, v1); + + double u_span = u1 - u0; + double v_span = v1 - v0; + + // 用步距估算扫描线数 + int num_passes = std::max(2, static_cast(v_span / std::max(step_over, 0.01)) + 1); + + for (int p = 0; p < num_passes; ++p) { + double v = v0 + v_span * static_cast(p) / static_cast(num_passes - 1); + for (int s = 0; s < samples_per_pass; ++s) { + double u = u0 + u_span * static_cast(s) / static_cast(samples_per_pass - 1); + Point3D pt = surface.evaluate(u, v); + Vector3D n = surface.normal(u, v); + result.emplace_back(pt, n); + } + } + return result; +} + +/// 从欧拉角计算刀轴方向 (ZYX 顺序，绕 X 为 tilt，绕 Y 为 lead) +Vector3D axis_from_angles(double tilt_deg, double lead_deg) { + double tx = deg2rad(tilt_deg); + double ty = deg2rad(lead_deg); + // 起始为 Z 轴 (0,0,1) + Vector3D axis(0, 0, 1); + // 绕 X 旋转 (tilt) + double ct = std::cos(tx), st = std::sin(tx); + double y1 = axis.y() * ct - axis.z() * st; + double z1 = axis.y() * st + axis.z() * ct; + // 绕 Y 旋转 (lead) + double cl = std::cos(ty), sl = std::sin(ty); + double x2 = z1 * sl + axis.x() * cl; + double z2 = z1 * cl - axis.x() * sl; + return Vector3D(x2, y1, z2).normalized(); +} + +/// 刀轴平滑插值（球面线性插值 Slerp） +Vector3D slerp_axis(const Vector3D& a, const Vector3D& b, double t) { + double dot = clamp(a.dot(b), -1.0, 1.0); + double theta = std::acos(dot); + if (theta < 1e-8) return a; + double st = std::sin(theta); + double w0 = std::sin((1.0 - t) * theta) / st; + double w1 = std::sin(t * theta) / st; + return (w0 * a + w1 * b).normalized(); +} + +/// 从向量计算旋转角（绕 X 的 tilt 和绕 Y 的 lead） +std::pair angles_from_axis(const Vector3D& axis) { + Vector3D n = axis.normalized(); + // tilt = asin(-n.y()) 绕 X 的旋转角 + double tilt = std::asin(clamp(-n.y(), -1.0, 1.0)); + // lead = atan2(n.x(), n.z()) 绕 Y 的旋转角 + double lead = std::atan2(n.x(), n.z()); + return {rad2deg(tilt), rad2deg(lead)}; +} + +} // anonymous namespace + +// ═══════════════════════════════════════════════════════════════════════════ +// MachineConfig +// ═══════════════════════════════════════════════════════════════════════════ + +bool MachineConfig::in_range(double a, double b, double c) const { + return a >= a_min && a <= a_max && + b >= b_min && b <= b_max && + c >= c_min && c <= c_max; +} + +std::optional> +MachineConfig::axis_to_angles(const Vector3D& tool_axis) const { + Vector3D n = tool_axis.normalized(); + + if (axis_config == MachineAxisConfig::AC) { + // A 绕 X, C 绕 Z + // tool_axis = Rz(C) * Rx(A) * (0,0,1) + // A = acos(n.z), C = atan2(n.y, n.x) — 简化 + double a = rad2deg(std::acos(clamp(n.z(), -1.0, 1.0))); + double c = rad2deg(std::atan2(n.y(), n.x())); + // 修正：更精确的 AC 分解 + if (std::abs(n.z()) > 0.9999) { + a = (n.z() > 0) ? 0.0 : 180.0; + c = 0.0; + } else { + a = rad2deg(std::acos(clamp(n.z(), -1.0, 1.0))); + c = rad2deg(std::atan2(n.y(), -n.x())); + } + c = normalize_angle(c); + if (in_range(a, 9999.0, c)) return std::make_tuple(a, 0.0, c); + // 尝试备选解 + double a2 = -a; + if (in_range(a2, 9999.0, normalize_angle(c + 180.0))) + return std::make_tuple(a2, 0.0, normalize_angle(c + 180.0)); + return std::nullopt; + } + + if (axis_config == MachineAxisConfig::BC) { + // B 绕 Y, C 绕 Z + if (std::abs(n.y()) > 0.9999) { + return std::make_tuple(0.0, (n.y() > 0) ? 90.0 : -90.0, 0.0); + } + double b = rad2deg(std::asin(clamp(-n.x(), -1.0, 1.0))); + double c = rad2deg(std::atan2(n.x(), n.z())); + c = normalize_angle(c); + if (in_range(9999.0, b, c)) return std::make_tuple(0.0, b, c); + return std::nullopt; + } + + if (axis_config == MachineAxisConfig::AB) { + // A 绕 X, B 绕 Y + if (std::abs(n.z()) > 0.9999) { + return std::make_tuple(0.0, 0.0, 0.0); + } + double a = rad2deg(std::asin(clamp(-n.y(), -1.0, 1.0))); + double b = rad2deg(std::asin(clamp(n.x(), -1.0, 1.0))); + if (in_range(a, b, 9999.0)) return std::make_tuple(a, b, 0.0); + return std::nullopt; + } + + return std::nullopt; +} + +// ═══════════════════════════════════════════════════════════════════════════ +// swarf_machining — 侧刃加工 +// ═══════════════════════════════════════════════════════════════════════════ + +AxisToolpath5 swarf_machining( + const curves::NurbsSurface& surface, + const Tool& tool, + const SwarfParams& params) { + + AxisToolpath5 tp; + tp.name = "Swarf Machining"; + tp.safe_height = params.safe_height; + + // 沿曲面参数线生成切触点 + auto contacts = generate_contact_points(surface, params.step_over, 30); + if (contacts.empty()) return tp; + + double tilt_rad = deg2rad(params.tilt_angle); + + for (const auto& [pt, normal] : contacts) { + // 刀尖位置：沿法向偏移刀具半径以侧刃贴合 + double radius = tool.diameter * 0.5; + Point3D pos = pt; + + // 刀轴沿曲面法线偏转 + // 计算切削方向（近似用 u 向切线） + Vector3D cut_dir(1, 0, 0); // 默认 X 方向 + + // 刀轴 = normal 绕 cut_dir 旋转 tilt_angle + Vector3D nnorm = normal.normalized(); + Vector3D cnorm = cut_dir.normalized(); + + // Rodrigues 旋转公式 + double ct = std::cos(tilt_rad); + double st = std::sin(tilt_rad); + Vector3D axis = ct * nnorm + st * (cnorm.cross(nnorm)) + + (1.0 - ct) * (cnorm.dot(nnorm)) * cnorm; + axis.normalize(); + + tp.points.emplace_back(pos, axis, params.feed_rate); + } + + return tp; +} + +// ═══════════════════════════════════════════════════════════════════════════ +// multi_axis_roughing — 多轴粗加工 (3+2 定位) +// ═══════════════════════════════════════════════════════════════════════════ + +AxisToolpath5 multi_axis_roughing( + const brep::BrepModel& body, + const Tool& tool, + const MultiAxisRoughingParams& params, + const MachineConfig& config) { + + AxisToolpath5 tp; + tp.name = "Multi-Axis Roughing"; + tp.safe_height = params.safe_height; + + // 3+2 定位：计算固定刀轴方向 + Vector3D fixed_axis = axis_from_angles(params.tilt_angle, params.lead_angle); + + // 验证刀轴方向在机床行程内 + auto angles = config.axis_to_angles(fixed_axis); + if (!angles.has_value()) { + // 回退到 Z 轴 + fixed_axis = Vector3D::UnitZ(); + } + + // 获取工件的 mesh 并计算 Z 范围 + auto mesh = body.to_mesh(0.2); + if (mesh.num_faces() == 0) return tp; + + auto [z_min, z_max] = mesh_z_range(mesh); + + // 分层粗加工 + double effective_step = params.step_down; + if (effective_step <= 0.0) effective_step = (z_max - z_min) * 0.1; + if (effective_step <= 1e-9) effective_step = 1.0; + + double current_z = z_max; + // 偏移刀具半径以考虑刀轴倾斜 + double radius = tool.diameter * 0.5; + double xy_offset = radius * std::sin(deg2rad(params.tilt_angle)); + + // 快速移动到安全高度 + tp.points.emplace_back( + Point3D(0, 0, params.safe_height), fixed_axis, 5000.0); + + while (current_z > z_min + 1e-9) { + double slice_z = std::max(z_min, current_z - effective_step); + + // 在当前 Z 层生成环绕刀路（XY 平面轮廓扫描） + // 简化：生成矩形包围盒的螺旋刀路 + AABB3D bbox = body.bounds(); + double x_min = bbox.min().x() - xy_offset - params.stock_to_leave; + double x_max = bbox.max().x() + xy_offset + params.stock_to_leave; + double y_min = bbox.min().y() - xy_offset - params.stock_to_leave; + double y_max = bbox.max().y() + xy_offset + params.stock_to_leave; + + double step = params.step_over; + if (step <= 0.0) step = (x_max - x_min) * 0.1; + + for (double x = x_min; x <= x_max; x += step) { + tp.points.emplace_back( + Point3D(x, y_min, slice_z), fixed_axis, params.feed_rate); + tp.points.emplace_back( + Point3D(x, y_max, slice_z), fixed_axis, params.feed_rate); + } + + current_z = slice_z; + } + + // 回到安全高度 + tp.points.emplace_back( + Point3D(0, 0, params.safe_height), fixed_axis, 5000.0); + + return tp; +} + +// ═══════════════════════════════════════════════════════════════════════════ +// multi_axis_finishing — 多轴精加工 +// ═══════════════════════════════════════════════════════════════════════════ + +AxisToolpath5 multi_axis_finishing( + const curves::NurbsSurface& surface, + const brep::BrepModel& body, + const Tool& tool, + const MultiAxisFinishingParams& params, + const MachineConfig& config) { + + AxisToolpath5 tp; + tp.name = "Multi-Axis Finishing"; + tp.safe_height = params.safe_height; + + // 沿曲面参数线生成切触点 + auto contacts = generate_contact_points(surface, params.step_over, 25); + if (contacts.empty()) return tp; + + Vector3D prev_axis(0, 0, 1); + bool first = true; + + for (size_t i = 0; i < contacts.size(); ++i) { + const auto& [pt, normal] = contacts[i]; + Vector3D nnorm = normal.normalized(); + + Vector3D tool_axis; + switch (params.axis_strategy) { + case ToolAxisStrategy::NormalToSurface: + tool_axis = nnorm; + break; + case ToolAxisStrategy::FixedAngle: + tool_axis = axis_from_angles(params.max_tilt, params.lead_angle); + break; + case ToolAxisStrategy::Tilted: { + // 侧倾：绕切削方向旋转 + Vector3D cut_dir(1, 0, 0); + if (i > 0) { + cut_dir = (pt - std::get<0>(contacts[i-1])).normalized(); + } + double tr = deg2rad(params.max_tilt); + Vector3D cn = cut_dir.normalized(); + tool_axis = std::cos(tr) * nnorm + + std::sin(tr) * (cn.cross(nnorm)) + + (1.0 - std::cos(tr)) * (cn.dot(nnorm)) * cn; + break; + } + case ToolAxisStrategy::LeadLag: { + // 导前/滞后 + Vector3D cut_dir(1, 0, 0); + if (i > 0) { + cut_dir = (pt - std::get<0>(contacts[i-1])).normalized(); + } + double lr = deg2rad(params.lead_angle); + Vector3D cn = cut_dir.normalized(); + Vector3D perp = nnorm.cross(cn); + tool_axis = std::cos(lr) * nnorm + + std::sin(lr) * (perp.cross(nnorm)); + break; + } + } + + tool_axis.normalize(); + + // 刀轴平滑变化 (slerp 插值) + if (!first) { + if (prev_axis.dot(tool_axis) < 0.999) { + tool_axis = slerp_axis(prev_axis, tool_axis, 0.7); + } + } + prev_axis = tool_axis; + first = false; + + // 碰撞检测 + if (params.collision_check) { + bool collision = check_collision(pt, tool_axis, body, tool, + params.shank_clearance, + params.holder_clearance); + if (collision) { + // 碰撞时调整为更保守的轴方向 + tool_axis = slerp_axis(tool_axis, Vector3D::UnitZ(), 0.3); + } + } + + tp.points.emplace_back(pt, tool_axis, params.feed_rate); + } + + return tp; +} + +// ═══════════════════════════════════════════════════════════════════════════ +// inverse_kinematics — 机床运动学逆解 +// ═══════════════════════════════════════════════════════════════════════════ + +AxisToolpath5 inverse_kinematics( + const AxisToolpath5& toolpath, + const MachineConfig& config) { + + AxisToolpath5 result; + result.name = toolpath.name + " (Machine Coords)"; + + for (const auto& cl : toolpath.points) { + // 解析刀轴矢量为旋转角 + auto angles = config.axis_to_angles(cl.tool_axis); + if (!angles.has_value()) { + // 不可达 — 跳过或使用最近可达方向 + continue; + } + + auto [a, b, c] = angles.value(); + + // 根据机床类型计算线性轴补偿 + Point3D pos = cl.position; + if (config.machine_type == MachineType::HeadHead) { + // 双摆头：刀具绕旋转中心旋转后偏移 + // 简化：用 RR^T * pivot 计算位置 + double ra = deg2rad(a); + double rb = deg2rad(b); + double rc = deg2rad(c); + + // 旋转中心在 pivot_offset + Vector3D pivot = config.pivot_offset; + Vector3D tool_tip(0, 0, -config.tool_length); + + // 近似：直接使用原始位置 + // 完整运动学：P_machine = P_work - R(a,b,c) * (pivot + tool_tip) + } + + if (config.machine_type == MachineType::TableTable) { + // 双转台：工件绕旋转中心旋转 + // P_table = R_inv * (P_work - pivot) + pivot + double ra = deg2rad(a); + double rb = deg2rad(b); + double rc = deg2rad(c); + + Vector3D shifted = cl.position - config.pivot_offset; + + // 逆旋转回机床坐标 + // 绕 Z 反向 + double cr = std::cos(-rc), sr = std::sin(-rc); + double x1 = shifted.x() * cr - shifted.y() * sr; + double y1 = shifted.x() * sr + shifted.y() * cr; + + // 绕 X 反向 + double ca = std::cos(-ra), sa = std::sin(-ra); + double y2 = y1 * ca - shifted.z() * sa; + double z2 = y1 * sa + shifted.z() * ca; + + pos = Point3D(x1, y2, z2) + config.pivot_offset; + } + + // 存储机床坐标：position = X/Y/Z, tool_axis 中存 A/B/C 度数 + CLData machine_cl(pos, Vector3D(a, b, c), cl.feed_rate); + result.points.push_back(machine_cl); + } + + return result; +} + +// ═══════════════════════════════════════════════════════════════════════════ +// check_collision — 刀柄/刀杆碰撞检测 +// ═══════════════════════════════════════════════════════════════════════════ + +bool check_collision( + const Point3D& position, + const Vector3D& tool_axis, + const brep::BrepModel& body, + const Tool& tool, + double shank_clearance, + double holder_clearance) { + + // 将工件 tessellate 为网格用于快速检测 + auto mesh = body.to_mesh(0.1); + if (mesh.num_faces() == 0) return false; + + Vector3D axis = tool_axis.normalized(); + double radius = tool.diameter * 0.5; + + // 检查刀柄区域 (从刃长结束到总长) + double shank_start = tool.length; + double shank_end = tool.overall_length; + + // 沿刀轴采样多个点检测间隙 + int samples = 8; + for (int i = 0; i <= samples; ++i) { + double t = static_cast(i) / static_cast(samples); + + // 刀柄圆柱采样圆上的点 + double z_offset = shank_start + t * (shank_end - shank_start); + double check_radius = (t < 0.5) + ? radius + shank_clearance // 近端：刀柄直径+间隙 + : tool.shank_diameter * 0.5 + holder_clearance; // 远端：刀杆 + + // 在垂直于刀轴的圆上采样 + Vector3D perp_x, perp_y; + if (std::abs(axis.z()) < 0.999) { + perp_x = axis.cross(Vector3D::UnitZ()).normalized(); + } else { + perp_x = axis.cross(Vector3D::UnitX()).normalized(); + } + perp_y = axis.cross(perp_x).normalized(); + + int angular_samples = 4; + for (int a = 0; a < angular_samples; ++a) { + double angle = 2.0 * M_PI * static_cast(a) / static_cast(angular_samples); + Vector3D offset = std::cos(angle) * perp_x + std::sin(angle) * perp_y; + Point3D sample_pt = position + z_offset * axis + check_radius * offset; + + // 对工件网格做射线检测 + double t_hit; + if (ray_mesh_intersect(sample_pt, -axis, mesh, z_offset, t_hit)) { + return true; // 碰撞 + } + } + } + + return false; +} + +// ═══════════════════════════════════════════════════════════════════════════ +// optimize_tool_axis — 刀轴方向优化 +// ═══════════════════════════════════════════════════════════════════════════ + +Vector3D optimize_tool_axis( + const std::vector& candidates, + const MachineConfig& config) { + + if (candidates.empty()) return Vector3D::UnitZ(); + + double best_score = -std::numeric_limits::max(); + Vector3D best_axis = Vector3D::UnitZ(); + + for (const auto& axis : candidates) { + Vector3D n = axis.normalized(); + + // 检查机床可达性 + auto angles = config.axis_to_angles(n); + if (!angles.has_value()) continue; + + auto [a, b, c] = angles.value(); + + // 评分：越接近 Z 轴分越高，同时旋转轴行程越小分越高 + double z_score = n.dot(Vector3D::UnitZ()); // [-1, 1], 接近 1 优先 + double angle_score = 1.0 - (std::abs(a) + std::abs(b)) / 180.0; + double score = z_score * 0.6 + angle_score * 0.4; + + // 避免刀尖切削（刀轴接近水平时扣分） + if (std::abs(n.z()) < 0.1) score -= 2.0; + + if (score > best_score) { + best_score = score; + best_axis = n; + } + } + + return best_axis; +} + +} // namespace vde::core diff --git a/src/mesh/fea_mesh.cpp b/src/mesh/fea_mesh.cpp new file mode 100644 index 0000000..1e74946 --- /dev/null +++ b/src/mesh/fea_mesh.cpp @@ -0,0 +1,938 @@ +#include "vde/mesh/fea_mesh.h" +#include "vde/mesh/mesh_quality.h" +#include "vde/brep/modeling.h" +#include "vde/core/aabb.h" +#include +#include +#include +#include +#include +#include +#include + +namespace vde::mesh { + +using core::AABB3D; + +// ═══════════════════════════════════════════════════════════ +// 内部辅助：从B-Rep tessellate 获取表面网格 +// ═══════════════════════════════════════════════════════════ + +namespace { + +/** 从 BrepModel 获取 tessellated 表面三角网格 */ +HalfedgeMesh tessellate_body(const brep::BrepModel& body, double deflection = 0.01) { + return body.to_mesh(deflection); +} + +/** 计算表面法向（用于边界层拉伸方向） */ +Vector3D compute_surface_normal(const HalfedgeMesh& mesh, int fi) { + return mesh.face_normal(fi); +} + +/** 计算实体包围盒 */ +AABB3D compute_bounds(const HalfedgeMesh& mesh) { + if (mesh.num_vertices() == 0) return AABB3D{}; + AABB3D bbox; + for (size_t i = 0; i < mesh.num_vertices(); ++i) { + bbox.expand(mesh.vertex(i)); + } + return bbox; +} + +/** 生成内部采样点（用于Delaunay四面体化） */ +std::vector generate_internal_points(const HalfedgeMesh& surface, + double target_size) { + auto bb = compute_bounds(surface); + Vector3D ext = {bb.max().x() - bb.min().x(), + bb.max().y() - bb.min().y(), + bb.max().z() - bb.min().z()}; + + if (target_size <= 0.0) { + target_size = std::max({ext.x(), ext.y(), ext.z()}) * 0.1; + } + + int nx = std::max(2, static_cast(ext.x() / target_size)); + int ny = std::max(2, static_cast(ext.y() / target_size)); + int nz = std::max(2, static_cast(ext.z() / target_size)); + + std::vector internal; + for (int i = 0; i < nx; ++i) { + for (int j = 0; j < ny; ++j) { + for (int k = 0; k < nz; ++k) { + double x = bb.min().x() + ext.x() * (i + 0.5) / nx; + double y = bb.min().y() + ext.y() * (j + 0.5) / ny; + double z = bb.min().z() + ext.z() * (k + 0.5) / nz; + internal.push_back(Point3D{x, y, z}); + } + } + } + return internal; +} + +// ═══════════════════════════════════════════════════════ +// 四面体网格质量优化 +// ═══════════════════════════════════════════════════════ + +/** Laplacian 光顺内部节点（不移动表面顶点） */ +void laplacian_smooth_tets(std::vector& verts, + const std::vector>& tets, + int surface_vertex_count, + int iterations) { + // 构建邻接表 + std::vector> neighbors(verts.size()); + for (auto& t : tets) { + for (int i = 0; i < 4; ++i) { + for (int j = i+1; j < 4; ++j) { + neighbors[t[i]].push_back(t[j]); + neighbors[t[j]].push_back(t[i]); + } + } + } + + for (int iter = 0; iter < iterations; ++iter) { + std::vector new_verts = verts; + for (size_t vi = static_cast(surface_vertex_count); vi < verts.size(); ++vi) { + auto& nbrs = neighbors[vi]; + if (nbrs.empty()) continue; + // 去重 + std::sort(nbrs.begin(), nbrs.end()); + nbrs.erase(std::unique(nbrs.begin(), nbrs.end()), nbrs.end()); + + Vector3D centroid{0,0,0}; + for (int n : nbrs) { + auto& np = verts[n]; + centroid = centroid + (np - Point3D{0,0,0}); + // centroid += np coords + centroid = {centroid.x() + np.x(), centroid.y() + np.y(), centroid.z() + np.z()}; + } + // Actually simpler: + } + } + + // Simpler implementation: + for (int iter = 0; iter < iterations; ++iter) { + std::vector new_verts = verts; + for (size_t vi = static_cast(surface_vertex_count); vi < verts.size(); ++vi) { + if (neighbors[vi].empty()) continue; + Vector3D sum{0,0,0}; + int count = 0; + for (int n : neighbors[vi]) { + sum = {sum.x() + verts[n].x(), + sum.y() + verts[n].y(), + sum.z() + verts[n].z()}; + ++count; + } + if (count > 0) { + new_verts[vi] = Point3D{sum.x()/count, sum.y()/count, sum.z()/count}; + } + } + verts.swap(new_verts); + } +} + +} // anonymous namespace + +// ═══════════════════════════════════════════════════════════ +// tetrahedral_mesh +// ═══════════════════════════════════════════════════════════ + +FEAMesh tetrahedral_mesh(const brep::BrepModel& body, + const TetMeshParams& params) { + // Step 1: Tessellate B-Rep surface + HalfedgeMesh surface = tessellate_body(body); + + // Step 2: Collect surface vertices + generate internal points + std::vector all_points; + int surface_vertex_count = static_cast(surface.num_vertices()); + for (size_t i = 0; i < surface.num_vertices(); ++i) { + all_points.push_back(surface.vertex(i)); + } + + double target_size = params.max_size; + if (target_size <= 0.0) { + auto bb = compute_bounds(surface); + Vector3D ext = {bb.max().x() - bb.min().x(), + bb.max().y() - bb.min().y(), + bb.max().z() - bb.min().z()}; + target_size = std::max({ext.x(), ext.y(), ext.z()}) * 0.05; + } + auto internal_pts = generate_internal_points(surface, target_size); + all_points.insert(all_points.end(), internal_pts.begin(), internal_pts.end()); + + // Step 3: Delaunay 3D 四面体化 + auto tet_mesh = delaunay_3d(all_points); + + // Step 4: 质量优化（Laplacian光顺内部节点） + if (params.quality_iterations > 0) { + laplacian_smooth_tets(tet_mesh.vertices, tet_mesh.tetrahedra, + surface_vertex_count, params.quality_iterations); + } + + // Step 5: 构建 FEAMesh + FEAMesh result; + result.element_type = FEAElementType::Tet4; + result.vertices = std::move(tet_mesh.vertices); + for (auto& t : tet_mesh.tetrahedra) { + result.elements.push_back({t[0], t[1], t[2], t[3]}); + } + result.material_ids.resize(result.elements.size(), 0); + + // Extract boundary faces (triangles on convex hull) + // For each face of each tet, if it appears only once, it's a boundary face + // Use string keys instead of array (C++17 compatible) + auto face_key = [](int a, int b, int c) -> std::string { + int mn = std::min({a,b,c}), mx = std::max({a,b,c}); + int md = a+b+c - mn - mx; + return std::to_string(mn) + "," + std::to_string(md) + "," + std::to_string(mx); + }; + + std::map face_count; + std::map> face_data; + + for (auto& t : tet_mesh.tetrahedra) { + // 4 faces of a tetrahedron + std::array, 4> tfaces = {{ + {t[0], t[1], t[2]}, {t[0], t[3], t[1]}, + {t[0], t[2], t[3]}, {t[1], t[3], t[2]} + }}; + for (auto& f : tfaces) { + auto key = face_key(f[0], f[1], f[2]); + face_count[key]++; + face_data[key] = f; + } + } + for (auto& kv : face_count) { + if (kv.second == 1) { // boundary face + auto& f = face_data[kv.first]; + std::vector bf; + bf.push_back(f[0]); bf.push_back(f[1]); bf.push_back(f[2]); + result.boundary_faces.push_back(bf); + } + } + result.boundary_face_tags.resize(result.boundary_faces.size(), 0); + + return result; +} + +// ═══════════════════════════════════════════════════════════ +// boundary_layer_mesh +// ═══════════════════════════════════════════════════════════ + +FEAMesh boundary_layer_mesh(const brep::BrepModel& body, + const BLPParams& params) { + // Step 1: Tessellate surface + HalfedgeMesh surface = tessellate_body(body); + + FEAMesh result; + result.element_type = FEAElementType::Wedge6; + + // Step 2: 计算每层高度 + int layers = params.num_layers; + std::vector heights(layers); + double h = params.first_cell_height; + double cumulative = 0; + for (int i = 0; i < layers; ++i) { + heights[i] = params.first_cell_height; + cumulative += h; + h *= params.growth_rate; + } + + // Step 3: 对每个表面顶点，沿法向生成 prism 层 + // vertex i of surface → prisms: base at original, top at base + normal*dist + size_t nv = surface.num_vertices(); + size_t nf = surface.num_faces(); + + // 计算顶点法向（面积加权平均） + std::vector vnormals(nv, Vector3D{0,0,0}); + for (size_t fi = 0; fi < nf; ++fi) { + Vector3D fn = surface.face_normal(static_cast(fi)); + auto fv = surface.face_vertices(static_cast(fi)); + for (int vi : fv) { + vnormals[vi] = {vnormals[vi].x() + fn.x(), + vnormals[vi].y() + fn.y(), + vnormals[vi].z() + fn.z()}; + } + } + for (auto& n : vnormals) { + double len = std::sqrt(n.x()*n.x() + n.y()*n.y() + n.z()*n.z()); + if (len > 1e-12) { + n = {n.x()/len, n.y()/len, n.z()/len}; + } + } + + // 生成顶点: 原始表面顶点 + 每层偏移顶点 + result.vertices.reserve(nv * (layers + 1)); + for (size_t vi = 0; vi < nv; ++vi) { + result.vertices.push_back(surface.vertex(vi)); + } + for (int layer = 0; layer < layers; ++layer) { + double dist = 0; + for (int l = 0; l <= layer; ++l) dist += heights[l]; + for (size_t vi = 0; vi < nv; ++vi) { + auto& vn = vnormals[vi]; + auto& sp = surface.vertex(vi); + result.vertices.push_back(Point3D{ + sp.x() + vn.x() * dist, + sp.y() + vn.y() * dist, + sp.z() + vn.z() * dist + }); + } + } + + // 生成棱柱单元 (Wedge6) + // 每个表面三角面 → layers 个 prism + // Wedge6 顺序: 底面(v0,v1,v2), 顶面(v3,v4,v5)，顶点对应: v3 在 v0 上方 + for (size_t fi = 0; fi < nf; ++fi) { + auto fv = surface.face_vertices(static_cast(fi)); + if (fv.size() != 3) continue; + int v0 = fv[0], v1 = fv[1], v2 = fv[2]; + + for (int layer = 0; layer < layers; ++layer) { + int base_offset = static_cast(nv * layer); + int top_offset = static_cast(nv * (layer + 1)); + // Wedge6: 底面 v0,v1,v2 → 顶面 v0',v1',v2' + result.elements.push_back({ + v0 + base_offset, v1 + base_offset, v2 + base_offset, // bottom + v0 + top_offset, v1 + top_offset, v2 + top_offset // top + }); + } + } + result.material_ids.resize(result.elements.size(), 0); + + // 提取边界（顶面三角） + int top_offset = static_cast(nv * layers); + for (size_t fi = 0; fi < nf; ++fi) { + auto fv = surface.face_vertices(static_cast(fi)); + if (fv.size() != 3) continue; + result.boundary_faces.push_back({ + fv[0] + top_offset, fv[1] + top_offset, fv[2] + top_offset + }); + } + result.boundary_face_tags.resize(result.boundary_faces.size(), 1); + + return result; +} + +// ═══════════════════════════════════════════════════════════ +// hexahedral_mesh +// ═══════════════════════════════════════════════════════════ + +FEAMesh hexahedral_mesh(const brep::BrepModel& body, + const HexMeshParams& params) { + // Step 1: Tessellate surface + HalfedgeMesh surface = tessellate_body(body); + + // Step 2: Sweep 方向 — 简化实现：沿 Z 轴 sweep + // 实际应用中应根据几何自动检测 sweep 方向 + auto bb = compute_bounds(surface); + double sweep_dist = params.sweep_distance; + if (sweep_dist <= 0.0) { + sweep_dist = bb.max().z() - bb.min().z(); + if (sweep_dist <= 0.0) sweep_dist = 1.0; + } + + int layers = std::max(1, params.sweep_layers); + double layer_height = sweep_dist / layers; + + FEAMesh result; + result.element_type = FEAElementType::Hex8; + + size_t nv = surface.num_vertices(); + size_t nf = surface.num_faces(); + + // 分离底面和顶面顶点（简化：Z ≈ z_min 为底面，Z ≈ z_max 为顶面） + std::vector bottom_verts, top_verts; + for (size_t vi = 0; vi < nv; ++vi) { + double z = surface.vertex(vi).z(); + if (std::abs(z - bb.min().z()) < 1e-6) { + bottom_verts.push_back(vi); + } else if (std::abs(z - bb.max().z()) < 1e-6) { + top_verts.push_back(vi); + } + } + + // 如果没有明显分层，使用所有顶点作为底面 + if (bottom_verts.empty()) { + for (size_t vi = 0; vi < nv; ++vi) bottom_verts.push_back(vi); + } + + // 生成顶点层 + result.vertices.reserve(nv * (layers + 1)); + for (size_t vi = 0; vi < nv; ++vi) { + result.vertices.push_back(Point3D{ + surface.vertex(vi).x(), + surface.vertex(vi).y(), + bb.min().z() + }); + } + for (int layer = 1; layer <= layers; ++layer) { + for (size_t vi = 0; vi < nv; ++vi) { + result.vertices.push_back(Point3D{ + surface.vertex(vi).x(), + surface.vertex(vi).y(), + bb.min().z() + layer * layer_height + }); + } + } + + // 生成六面体单元: 对每个表面四边形边界面生成 sweep hex + // 简化: 使用三角形 split 成四边形 + for (int layer = 0; layer < layers; ++layer) { + int base_off = static_cast(nv * layer); + int top_off = static_cast(nv * (layer + 1)); + + for (size_t fi = 0; fi < nf; ++fi) { + auto fv = surface.face_vertices(static_cast(fi)); + if (fv.size() == 4) { + // 四边形面 → 直接 sweep + result.elements.push_back({ + fv[0]+base_off, fv[1]+base_off, fv[2]+base_off, fv[3]+base_off, + fv[0]+top_off, fv[1]+top_off, fv[2]+top_off, fv[3]+top_off + }); + } else if (fv.size() == 3 && params.use_quad_dominant) { + // 三角形 → 拆分为3个四边形（连接各边中点到中心） + // 简化：将三角形视为退化的四边形（v3 = v2） + result.elements.push_back({ + fv[0]+base_off, fv[1]+base_off, fv[2]+base_off, fv[2]+base_off, + fv[0]+top_off, fv[1]+top_off, fv[2]+top_off, fv[2]+top_off + }); + } + } + } + result.material_ids.resize(result.elements.size(), 0); + result.boundary_face_tags.clear(); + + return result; +} + +// ═══════════════════════════════════════════════════════════ +// adaptive_refinement +// ═══════════════════════════════════════════════════════════ + +FEAMesh adaptive_refinement(const FEAMesh& mesh, + const ErrorEstimator& estimator, + const RefinementParams& params) { + FEAMesh result = mesh; + + for (int iter = 0; iter < params.max_iterations; ++iter) { + // Step 1: 评估所有单元的误差 + std::vector errors(result.elements.size()); + double max_error = 0.0; + for (size_t ei = 0; ei < result.elements.size(); ++ei) { + errors[ei] = estimator(static_cast(ei), result); + max_error = std::max(max_error, errors[ei]); + } + if (max_error < params.error_threshold * 1e-3) break; + + // Step 2: 标记需要细化的单元 + double threshold = max_error * params.error_threshold; + std::vector refine_mark(result.elements.size(), false); + size_t refine_count = 0; + for (size_t ei = 0; ei < result.elements.size(); ++ei) { + if (errors[ei] > threshold) { + refine_mark[ei] = true; + ++refine_count; + } + } + if (refine_count == 0) break; + + // Step 3: 边中点插入细分（Tet4 → 子四面体） + // 使用边的中点将每个标记的四面体细分为4个小子四面体 + std::vector new_verts = result.vertices; + std::vector> new_elements; + + // 边中点缓存: (v_i, v_j) → new_vertex_index + std::map, int> edge_midpoints; + + auto get_edge_mid = [&](int a, int b) -> int { + if (a > b) std::swap(a, b); + auto key = std::make_pair(a, b); + auto it = edge_midpoints.find(key); + if (it != edge_midpoints.end()) return it->second; + int idx = static_cast(new_verts.size()); + Point3D pa = result.vertices[a], pb = result.vertices[b]; + new_verts.push_back(Point3D{ + (pa.x() + pb.x()) / 2.0, + (pa.y() + pb.y()) / 2.0, + (pa.z() + pb.z()) / 2.0 + }); + edge_midpoints[key] = idx; + return idx; + }; + + for (size_t ei = 0; ei < result.elements.size(); ++ei) { + if (!refine_mark[ei]) { + new_elements.push_back(result.elements[ei]); + continue; + } + + auto& e = result.elements[ei]; + if (e.size() < 4) { + new_elements.push_back(e); + continue; + } + + // Tet4 → 4 sub-tets (by splitting at centroid of edges) + // 简化: 使用6个边中点 + 原4个顶点 + int v0 = e[0], v1 = e[1], v2 = e[2], v3 = e[3]; + int m01 = get_edge_mid(v0, v1); + int m02 = get_edge_mid(v0, v2); + int m03 = get_edge_mid(v0, v3); + int m12 = get_edge_mid(v1, v2); + int m13 = get_edge_mid(v1, v3); + int m23 = get_edge_mid(v2, v3); + + // Subdivide into 4 corner tets + central octahedron (split into 4 tets) + // Corner tet at v0: {v0, m01, m02, m03} + new_elements.push_back({v0, m01, m02, m03}); + // Corner tet at v1: {v1, m01, m12, m13} + new_elements.push_back({v1, m01, m12, m13}); + // Corner tet at v2: {v2, m02, m12, m23} + new_elements.push_back({v2, m02, m12, m23}); + // Corner tet at v3: {v3, m03, m13, m23} + new_elements.push_back({v3, m03, m13, m23}); + // Remaining octahedron → 4 tets (choose the shorter diagonal) + new_elements.push_back({m01, m02, m12, m23}); + new_elements.push_back({m01, m02, m23, m03}); + new_elements.push_back({m01, m12, m23, m13}); + new_elements.push_back({m01, m02, m03, m13}); + // Actually this over-splits. Let me simplify: just 4 sub-tets per refinement + // Keep only the 4 corner tets for simplicity + } + + // Simplify: use only corner splits + new_elements.clear(); + edge_midpoints.clear(); + new_verts = result.vertices; + + for (size_t ei = 0; ei < result.elements.size(); ++ei) { + auto& e = result.elements[ei]; + if (!refine_mark[ei] || e.size() < 4) { + new_elements.push_back(e); + continue; + } + + int v0 = e[0], v1 = e[1], v2 = e[2], v3 = e[3]; + // 质量较高的简化细分：每个四面体 → 4个子四面体 + // 使用面重心 + 体心分裂 + // Face centroids + auto centroid3 = [&](int a, int b, int c) -> int { + if (a > b) std::swap(a,b); + if (b > c) std::swap(b,c); + if (a > b) std::swap(a,b); + auto key = std::make_pair(std::make_pair(a,b).first, c); + // use string key instead + auto key_str = std::to_string(a) + "," + std::to_string(b) + "," + std::to_string(c); + // Simple approach: just compute + Point3D pa = result.vertices[a], pb = result.vertices[b], pc = result.vertices[c]; + int idx = static_cast(new_verts.size()); + new_verts.push_back(Point3D{ + (pa.x()+pb.x()+pc.x())/3.0, + (pa.y()+pb.y()+pc.y())/3.0, + (pa.z()+pb.z()+pc.z())/3.0 + }); + return idx; + }; + + int cf0 = centroid3(v0, v1, v2); + int cf1 = centroid3(v0, v1, v3); + int cf2 = centroid3(v0, v2, v3); + int cf3 = centroid3(v1, v2, v3); + + // 4 sub-tets using face centroids + new_elements.push_back({v0, cf0, cf1, cf2}); + new_elements.push_back({v1, cf0, cf1, cf3}); + new_elements.push_back({v2, cf0, cf2, cf3}); + new_elements.push_back({v3, cf1, cf2, cf3}); + } + + result.vertices = std::move(new_verts); + result.elements = std::move(new_elements); + result.material_ids.resize(result.elements.size(), 0); + } + + return result; +} + +// ═══════════════════════════════════════════════════════════ +// 质量评估 — 单元素指标 +// ═══════════════════════════════════════════════════════════ + +double element_skewness(const std::vector& verts, FEAElementType type) { + // Skewness = (optimal_cell_size - actual_cell_size) / optimal_cell_size + // 对于 Tet: 基于等边四面体比较 + if (verts.size() < 4) return 1.0; + + auto tet_skewness = [](const Point3D& a, const Point3D& b, + const Point3D& c, const Point3D& d) -> double { + Vector3D e1 = b - a, e2 = c - a, e3 = d - a; + double vol = std::abs(e1.x()*(e2.y()*e3.z()-e2.z()*e3.y()) + -e1.y()*(e2.x()*e3.z()-e2.z()*e3.x()) + +e1.z()*(e2.x()*e3.y()-e2.y()*e3.x())) / 6.0; + double L1 = (b-a).norm(), L2 = (c-a).norm(), L3 = (d-a).norm(); + double L4 = (c-b).norm(), L5 = (d-b).norm(), L6 = (d-c).norm(); + double avg_L = (L1+L2+L3+L4+L5+L6)/6.0; + if (avg_L < 1e-12) return 1.0; + double opt_vol = avg_L*avg_L*avg_L / (6.0 * std::sqrt(2.0)); + double skew = (opt_vol - vol) / opt_vol; + return std::max(0.0, std::min(1.0, skew)); + }; + + if (type == FEAElementType::Tet4 && verts.size() >= 4) { + return tet_skewness(verts[0], verts[1], verts[2], verts[3]); + } + + auto tet_skewness_hex = [&tet_skewness](const std::vector& v) -> double { + if (v.size() < 8) return 1.0; + // Decompose hex into 5(or 6) tets and average + double s = 0; + s += tet_skewness(v[0],v[1],v[3],v[4]); + s += tet_skewness(v[1],v[2],v[3],v[6]); + s += tet_skewness(v[1],v[3],v[4],v[6]); + s += tet_skewness(v[1],v[4],v[5],v[6]); + s += tet_skewness(v[3],v[4],v[6],v[7]); + return s / 5.0; + }; + + if (type == FEAElementType::Hex8) return tet_skewness_hex(verts); + if (type == FEAElementType::Wedge6 && verts.size() >= 6) { + double s = 0; + s += tet_skewness(verts[0],verts[1],verts[2],verts[3]); + s += tet_skewness(verts[1],verts[2],verts[3],verts[4]); + s += tet_skewness(verts[2],verts[3],verts[4],verts[5]); + return s / 3.0; + } + return 1.0; +} + +double element_aspect_ratio(const std::vector& verts, FEAElementType type) { + if (verts.size() < 2) return 1.0; + + // Compute all edge lengths + std::vector edge_lens; + for (size_t i = 0; i < verts.size(); ++i) { + for (size_t j = i+1; j < verts.size(); ++j) { + edge_lens.push_back((verts[j] - verts[i]).norm()); + } + } + if (edge_lens.empty()) return 1.0; + double min_e = *std::min_element(edge_lens.begin(), edge_lens.end()); + double max_e = *std::max_element(edge_lens.begin(), edge_lens.end()); + if (min_e < 1e-12) return 1e6; + return max_e / min_e; +} + +double element_scaled_jacobian(const std::vector& verts, FEAElementType type) { + if (verts.size() < 4) return 0.0; + + auto tet_jac = [](const Point3D& a, const Point3D& b, + const Point3D& c, const Point3D& d) -> double { + Vector3D e1 = b - a, e2 = c - a, e3 = d - a; + double J = e1.x()*(e2.y()*e3.z()-e2.z()*e3.y()) + -e1.y()*(e2.x()*e3.z()-e2.z()*e3.x()) + +e1.z()*(e2.x()*e3.y()-e2.y()*e3.x()); + double norm = e1.norm() * e2.norm() * e3.norm(); + if (norm < 1e-12) return 0.0; + return std::abs(J) / (norm * std::sqrt(2.0)/3.0); + }; + + if (type == FEAElementType::Tet4 && verts.size() >= 4) { + return tet_jac(verts[0], verts[1], verts[2], verts[3]); + } + + if (type == FEAElementType::Hex8 && verts.size() >= 8) { + // 5-tet decomposition of hex (standard scheme, no degenerate sub-tets) + // Hex nodes: bottom{0,1,2,3} CCW, top{4,5,6,7} CCW + double min_j = 1e6; + auto& v = verts; + min_j = std::min(min_j, tet_jac(v[0],v[1],v[2],v[5])); + min_j = std::min(min_j, tet_jac(v[0],v[2],v[7],v[5])); + min_j = std::min(min_j, tet_jac(v[0],v[2],v[3],v[7])); + min_j = std::min(min_j, tet_jac(v[0],v[5],v[7],v[4])); + min_j = std::min(min_j, tet_jac(v[2],v[7],v[5],v[6])); + return std::max(0.0, std::min(1.0, min_j)); + } + + if (type == FEAElementType::Wedge6 && verts.size() >= 6) { + double j = 0; + j += tet_jac(verts[0],verts[1],verts[2],verts[3]); + j += tet_jac(verts[0],verts[2],verts[4],verts[5]); + j += tet_jac(verts[0],verts[3],verts[4],verts[5]); + return std::max(0.0, std::min(1.0, j/3.0)); + } + + return 0.0; +} + +double element_orthogonality(const std::vector& verts, FEAElementType type) { + if (verts.size() < 4) return 0.0; + + auto tri_ortho = [](const Point3D& a, const Point3D& b, const Point3D& c, + const Point3D& o) -> double { + // Face (a,b,c) normal vs vector from face centroid to opposite vertex o + Vector3D e1 = b - a, e2 = c - a; + Vector3D n = {e1.y()*e2.z() - e1.z()*e2.y(), + e1.z()*e2.x() - e1.x()*e2.z(), + e1.x()*e2.y() - e1.y()*e2.x()}; + double nl = std::sqrt(n.x()*n.x()+n.y()*n.y()+n.z()*n.z()); + Point3D centroid{(a.x()+b.x()+c.x())/3, (a.y()+b.y()+c.y())/3, (a.z()+b.z()+c.z())/3}; + Vector3D vec = {o.x()-centroid.x(), o.y()-centroid.y(), o.z()-centroid.z()}; + double vl = std::sqrt(vec.x()*vec.x()+vec.y()*vec.y()+vec.z()*vec.z()); + if (nl < 1e-12 || vl < 1e-12) return 0.0; + double dot = (n.x()*vec.x()+n.y()*vec.y()+n.z()*vec.z())/(nl*vl); + return std::abs(dot); + }; + + if (type == FEAElementType::Tet4 && verts.size() >= 4) { + double o = 0; + o += tri_ortho(verts[1],verts[2],verts[3],verts[0]); + o += tri_ortho(verts[0],verts[2],verts[3],verts[1]); + o += tri_ortho(verts[0],verts[1],verts[3],verts[2]); + o += tri_ortho(verts[0],verts[1],verts[2],verts[3]); + return o / 4.0; + } + + if (type == FEAElementType::Hex8 && verts.size() >= 8) { + double o = 0; + o += tri_ortho(verts[4],verts[5],verts[6],verts[0]); // face 456 vs centroid to v0 + o += tri_ortho(verts[0],verts[1],verts[2],verts[6]); // opposite face + o += tri_ortho(verts[0],verts[1],verts[5],verts[3]); + o += tri_ortho(verts[1],verts[2],verts[6],verts[4]); + o += tri_ortho(verts[2],verts[3],verts[7],verts[5]); + o += tri_ortho(verts[0],verts[3],verts[7],verts[1]); + return o / 6.0; + } + + if (type == FEAElementType::Wedge6 && verts.size() >= 6) { + // 3 quad faces + 2 tri faces + double o = 0; + o += tri_ortho(verts[0],verts[1],verts[2],verts[3]); // tri bottom vs v3 + o += tri_ortho(verts[3],verts[4],verts[5],verts[0]); // tri top vs v0 + return o / 2.0; + } + + return 0.0; +} + +// ═══════════════════════════════════════════════════════════ +// fea_quality_report +// ═══════════════════════════════════════════════════════════ + +FEAQualityReport fea_quality_report(const FEAMesh& mesh) { + FEAQualityReport r; + r.total_elements = mesh.elements.size(); + if (r.total_elements == 0) return r; + + std::vector skews, aspects, jacobs, orthos; + + for (size_t ei = 0; ei < mesh.elements.size(); ++ei) { + // Collect element vertices + std::vector everts; + for (int vi : mesh.elements[ei]) { + everts.push_back(mesh.vertices[vi]); + } + + double sj = element_scaled_jacobian(everts, mesh.element_type); + jacobs.push_back(sj); + + if (sj < 0.0) r.inverted_elements++; + if (sj < 0.01) r.degenerate_elements++; + + skews.push_back(element_skewness(everts, mesh.element_type)); + aspects.push_back(element_aspect_ratio(everts, mesh.element_type)); + orthos.push_back(element_orthogonality(everts, mesh.element_type)); + } + + // Statistics + auto stats = [](const std::vector& v, double& mn, double& mx, double& avg) { + if (v.empty()) { mn = mx = avg = 0; return; } + mn = *std::min_element(v.begin(), v.end()); + mx = *std::max_element(v.begin(), v.end()); + avg = std::accumulate(v.begin(), v.end(), 0.0) / v.size(); + }; + + stats(skews, r.min_skewness, r.max_skewness, r.avg_skewness); + stats(aspects, r.min_aspect_ratio, r.max_aspect_ratio, r.avg_aspect_ratio); + stats(jacobs, r.min_jacobian, r.max_jacobian, r.avg_jacobian); + stats(orthos, r.min_orthogonality, r.max_orthogonality, r.avg_orthogonality); + + // Grades (by Jacobian) + for (auto j : jacobs) { + if (j >= 0.9) r.grade_A++; + else if (j >= 0.7) r.grade_B++; + else if (j >= 0.5) r.grade_C++; + else if (j >= 0.3) r.grade_D++; + else r.grade_F++; + } + + return r; +} + +// ═══════════════════════════════════════════════════════════ +// CAE 格式导出 +// ═══════════════════════════════════════════════════════════ + +namespace { + +/** Abaqus 单元类型名称 */ +const char* abaqus_element_name(FEAElementType t) { + switch (t) { + case FEAElementType::Tet4: return "C3D4"; + case FEAElementType::Tet10: return "C3D10"; + case FEAElementType::Hex8: return "C3D8"; + case FEAElementType::Hex20: return "C3D20"; + case FEAElementType::Wedge6: return "C3D6"; + case FEAElementType::Wedge15:return "C3D15"; + } + return "C3D4"; +} + +/** Nastran 单元类型名称 */ +const char* nastran_element_name(FEAElementType t) { + switch (t) { + case FEAElementType::Tet4: return "CTETRA"; + case FEAElementType::Tet10: return "CTETRA"; + case FEAElementType::Hex8: return "CHEXA"; + case FEAElementType::Hex20: return "CHEXA"; + case FEAElementType::Wedge6: return "CPENTA"; + case FEAElementType::Wedge15:return "CPENTA"; + } + return "CTETRA"; +} + +} // anonymous namespace + +bool export_abaqus(const FEAMesh& mesh, const std::string& filepath) { + FILE* fp = std::fopen(filepath.c_str(), "w"); + if (!fp) return false; + + std::fprintf(fp, "*HEADING\n"); + std::fprintf(fp, "** Generated by ViewDesignEngine FEA Mesh Module\n"); + std::fprintf(fp, "*NODE\n"); + + // 节点 + for (size_t i = 0; i < mesh.vertices.size(); ++i) { + auto& v = mesh.vertices[i]; + std::fprintf(fp, "%zu, %.10f, %.10f, %.10f\n", i+1, v.x(), v.y(), v.z()); + } + + // 单元 + const char* elm_name = abaqus_element_name(mesh.element_type); + std::fprintf(fp, "*ELEMENT, TYPE=%s, ELSET=ALL\n", elm_name); + for (size_t ei = 0; ei < mesh.elements.size(); ++ei) { + std::fprintf(fp, "%zu", ei+1); + for (int vi : mesh.elements[ei]) { + std::fprintf(fp, ", %d", vi+1); + } + std::fprintf(fp, "\n"); + } + + // 材料 + std::fprintf(fp, "*MATERIAL, NAME=MAT1\n"); + std::fprintf(fp, "*ELASTIC\n"); + std::fprintf(fp, "210000.0, 0.3\n"); + std::fprintf(fp, "*SOLID SECTION, ELSET=ALL, MATERIAL=MAT1\n"); + + // 边界面（作为 surface） + if (!mesh.boundary_faces.empty()) { + std::fprintf(fp, "*SURFACE, NAME=BOUNDARY, TYPE=ELEMENT\n"); + for (size_t bi = 0; bi < mesh.boundary_faces.size(); ++bi) { + auto& bf = mesh.boundary_faces[bi]; + std::fprintf(fp, "%zu, S%d\n", bi+1, bi+1); + } + } + + std::fprintf(fp, "*END ASSEMBLY\n"); + std::fclose(fp); + return true; +} + +bool export_ansys(const FEAMesh& mesh, const std::string& filepath) { + FILE* fp = std::fopen(filepath.c_str(), "w"); + if (!fp) return false; + + std::fprintf(fp, "/COM, Generated by ViewDesignEngine FEA Mesh Module\n"); + std::fprintf(fp, "/PREP7\n"); + + // 节点: NBLOCK format + std::fprintf(fp, "NBLOCK,6,SOLID\n"); + std::fprintf(fp, "(3i9,6e21.13e3)\n"); + for (size_t i = 0; i < mesh.vertices.size(); ++i) { + auto& v = mesh.vertices[i]; + // Ansys NBLOCK: id, solid, 0, 0, 0, x, y, z + std::fprintf(fp, "%9zu%9d%9d%21.13e%21.13e%21.13e\n", + i+1, 0, 0, v.x(), v.y(), v.z()); + } + std::fprintf(fp, "-1\n"); + + // 单元: EBLOCK format + int npe = mesh.npe(); + const char* elm_name = abaqus_element_name(mesh.element_type); + std::fprintf(fp, "EBLOCK,19,SOLID,%zu\n", mesh.elements.size()); + std::fprintf(fp, "(19i9)\n"); + + // Ansys element type mapping + int ansys_type = 185; // SOLID185 for hex + if (mesh.element_type == FEAElementType::Tet4) ansys_type = 185; + else if (mesh.element_type == FEAElementType::Hex8) ansys_type = 185; + else if (mesh.element_type == FEAElementType::Wedge6) ansys_type = 185; + + for (size_t ei = 0; ei < mesh.elements.size(); ++ei) { + std::fprintf(fp, "%9d%9d%9d%9d", 1, ansys_type, 1, 0); + for (int vi : mesh.elements[ei]) { + std::fprintf(fp, "%9d", vi+1); + } + // Pad remaining nodes with last node + for (int i = npe; i < 20; ++i) { + std::fprintf(fp, "%9d", mesh.elements[ei].back()+1); + } + std::fprintf(fp, "\n"); + } + std::fprintf(fp, "-1\n"); + + std::fprintf(fp, "FINISH\n"); + std::fclose(fp); + return true; +} + +bool export_nastran(const FEAMesh& mesh, const std::string& filepath) { + FILE* fp = std::fopen(filepath.c_str(), "w"); + if (!fp) return false; + + std::fprintf(fp, "$ Generated by ViewDesignEngine FEA Mesh Module\n"); + std::fprintf(fp, "BEGIN BULK\n"); + + // GRID cards + for (size_t i = 0; i < mesh.vertices.size(); ++i) { + auto& v = mesh.vertices[i]; + std::fprintf(fp, "GRID %-8zu %.6f %.6f %.6f\n", + i+1, v.x(), v.y(), v.z()); + } + + // Element cards + const char* elm_name = nastran_element_name(mesh.element_type); + for (size_t ei = 0; ei < mesh.elements.size(); ++ei) { + std::fprintf(fp, "%-8s %-8zu 0 ", elm_name, ei+1); + auto& e = mesh.elements[ei]; + for (int vi : e) { + std::fprintf(fp, "%-8d ", vi+1); + } + // Pad for Hex8 → CHEXA needs 6+ faces + if (mesh.element_type == FEAElementType::Hex8 && e.size() == 8) { + std::fprintf(fp, "\n "); + // Second line for CHEXA: nodes 7-8 (already included above if 8 nodes) + // Actually Nastran CHEXA: nodes 1-8 on first line, nothing needed + } + std::fprintf(fp, "\n"); + } + + // Material + std::fprintf(fp, "MAT1 1 210000.0 0.3\n"); + + std::fprintf(fp, "ENDDATA\n"); + std::fclose(fp); + return true; +} + +} // namespace vde::mesh diff --git a/src/mesh/point_cloud.cpp b/src/mesh/point_cloud.cpp new file mode 100644 index 0000000..57a1022 --- /dev/null +++ b/src/mesh/point_cloud.cpp @@ -0,0 +1,731 @@ +#include "vde/mesh/point_cloud.h" +#include +#include +#include +#include +#include +#include +#include +#include + +namespace vde::mesh { + +// ─────────────────────────────────────────────────────────── +// bounds +// ─────────────────────────────────────────────────────────── + +AABB3D PointCloud::bounds() const { + AABB3D bbox; + for (const auto& p : vertices) bbox.expand(p); + return bbox; +} + +// ─────────────────────────────────────────────────────────── +// PLY I/O +// ─────────────────────────────────────────────────────────── + +PointCloud PointCloud::load_ply(const std::string& path) { + std::ifstream in(path, std::ios::binary); + if (!in) throw std::runtime_error("Cannot open PLY file: " + path); + + std::string line; + size_t num_vertices = 0; + bool ascii = true; + struct Prop { std::string name; bool is_uchar = false; }; + std::vector props; + int x_idx = -1, y_idx = -1, z_idx = -1; + int nx_idx = -1, ny_idx = -1, nz_idx = -1; + int r_idx = -1, g_idx = -1, b_idx = -1; + // Normal indices tracked but extracted in binary mode below; suppress unused warning + (void)nx_idx; (void)ny_idx; (void)nz_idx; + + // Read header + while (std::getline(in, line)) { + if (line.rfind("format ascii", 0) == 0) { ascii = true; } + else if (line.rfind("format binary_little_endian", 0) == 0) { ascii = false; } + else if (line.rfind("element vertex ", 0) == 0) { + num_vertices = std::stoul(line.substr(15)); + } + else if (line.rfind("property ", 0) == 0) { + std::istringstream iss(line.substr(9)); + std::string type, name; + iss >> type >> name; + if (type == "float" || type == "float32" || type == "double") { + props.push_back({name, false}); + } else if (type == "uchar" || type == "uint8") { + props.push_back({name, true}); + } else { + props.push_back({name, false}); + } + } + else if (line == "end_header") break; + } + + // Map properties + for (size_t i = 0; i < props.size(); ++i) { + if (props[i].name == "x") x_idx = static_cast(i); + if (props[i].name == "y") y_idx = static_cast(i); + if (props[i].name == "z") z_idx = static_cast(i); + if (props[i].name == "nx") nx_idx = static_cast(i); + if (props[i].name == "ny") ny_idx = static_cast(i); + if (props[i].name == "nz") nz_idx = static_cast(i); + if (props[i].name == "red" || props[i].name == "r" || props[i].name == "diffuse_red") r_idx = static_cast(i); + if (props[i].name == "green" || props[i].name == "g" || props[i].name == "diffuse_green") g_idx = static_cast(i); + if (props[i].name == "blue" || props[i].name == "b" || props[i].name == "diffuse_blue") b_idx = static_cast(i); + } + + PointCloud cloud; + cloud.vertices.reserve(num_vertices); + + if (ascii) { + for (size_t i = 0; i < num_vertices; ++i) { + std::getline(in, line); + std::istringstream iss(line); + std::vector vals(props.size()); + for (size_t j = 0; j < props.size(); ++j) iss >> vals[j]; + + if (x_idx >= 0 && y_idx >= 0 && z_idx >= 0) + cloud.vertices.emplace_back(vals[x_idx], vals[y_idx], vals[z_idx]); + } + } else { + // Binary little-endian + // Calculate per-vertex byte size + size_t vertex_size = 0; + for (const auto& p : props) vertex_size += p.is_uchar ? 1 : 4; + + std::vector buf(vertex_size * num_vertices); + in.read(buf.data(), static_cast(buf.size())); + + for (size_t i = 0; i < num_vertices; ++i) { + const char* base = buf.data() + i * vertex_size; + size_t off = 0; + float x = 0, y = 0, z = 0; + for (size_t j = 0; j < props.size(); ++j) { + if (props[j].is_uchar) { + uint8_t val = static_cast(base[off]); + off += 1; + if (static_cast(j) == r_idx) { /* store color later */ } + if (static_cast(j) == g_idx) { /* store color later */ } + if (static_cast(j) == b_idx) { /* store color later */ } + } else { + float val; + std::memcpy(&val, base + off, 4); + off += 4; + if (static_cast(j) == x_idx) x = val; + if (static_cast(j) == y_idx) y = val; + if (static_cast(j) == z_idx) z = val; + } + } + cloud.vertices.emplace_back(x, y, z); + } + } + + return cloud; +} + +bool PointCloud::save_ply(const std::string& path) const { + std::ofstream out(path); + if (!out) return false; + + out << "ply\nformat ascii 1.0\n"; + out << "element vertex " << vertices.size() << "\n"; + out << "property float x\nproperty float y\nproperty float z\n"; + if (has_normals()) { + out << "property float nx\nproperty float ny\nproperty float nz\n"; + } + if (has_colors()) { + out << "property uchar red\nproperty uchar green\nproperty uchar blue\n"; + } + out << "end_header\n"; + + for (size_t i = 0; i < vertices.size(); ++i) { + const auto& v = vertices[i]; + out << v.x() << " " << v.y() << " " << v.z(); + if (has_normals()) { + out << " " << normals[i].x() << " " << normals[i].y() << " " << normals[i].z(); + } + if (has_colors()) { + out << " " << static_cast(colors[i][0]) + << " " << static_cast(colors[i][1]) + << " " << static_cast(colors[i][2]); + } + out << "\n"; + } + return true; +} + +// ─────────────────────────────────────────────────────────── +// E57 (simplified XML+Binary parser) +// ─────────────────────────────────────────────────────────── + +PointCloud PointCloud::load_e57(const std::string& path) { + // E57 is a complex binary format combining XML metadata + binary sections. + // Full spec: ASTM E2807. Here we provide a simplified parser that + // extracts CartesianBounds + points3D sections. + // + // The E57 file is structured as: + // - Header (48 bytes): magic, version, XML offset, XML length + // - XML section: metadata describing data3D structures + // - Binary section: packed point data + // + // For production use, use libE57 or PCL. This simplified parser + // handles basic Leica/FARO scans. + + std::ifstream in(path, std::ios::binary); + if (!in) throw std::runtime_error("Cannot open E57 file: " + path); + + // Read header + char header[48]; + in.read(header, 48); + if (std::string(header, 8) != "ASTM-E57") { + throw std::runtime_error("Not a valid E57 file (bad magic): " + path); + } + + // Parse simplified XML to find points array info + uint64_t xml_offset = 0, xml_length = 0; + std::memcpy(&xml_offset, header + 24, 8); + std::memcpy(&xml_length, header + 32, 8); + + in.seekg(static_cast(xml_offset)); + std::string xml_data(xml_length, '\0'); + in.read(&xml_data[0], static_cast(xml_length)); + + // Find pointCount in XML + size_t point_count = 0; + { + auto pos = xml_data.find("pointCount"); + if (pos != std::string::npos) { + auto val_start = xml_data.find('>', pos); + auto val_end = xml_data.find('<', val_start); + if (val_start != std::string::npos && val_end != std::string::npos) + point_count = std::stoull(xml_data.substr(val_start + 1, val_end - val_start - 1)); + } + } + + if (point_count == 0) { + // Fallback: estimate from file size + point_count = (xml_offset - 48) / (3 * 8 + 4); // x,y,z doubles + intensity float + } + + // Find binary data offset (after XML) + uint64_t binary_offset = xml_offset + xml_length; + // binary data: pointCount records of {double x, double y, double z, float intensity, uint8_t r, uint8_t g, uint8_t b} + // For simplicity, read the 3 doubles (24 bytes) per point + in.seekg(static_cast(binary_offset)); + + PointCloud cloud; + cloud.vertices.reserve(point_count); + + const size_t record_size = 3 * 8; // x, y, z as doubles + for (size_t i = 0; i < point_count; ++i) { + double xyz[3]; + if (!in.read(reinterpret_cast(xyz), record_size)) break; + cloud.vertices.emplace_back(xyz[0], xyz[1], xyz[2]); + } + + return cloud; +} + +// ─────────────────────────────────────────────────────────── +// PTX (Leica scanner format) +// ─────────────────────────────────────────────────────────── + +PointCloud PointCloud::load_ptx(const std::string& path) { + // PTX format: text format, each scan section has: + // Line 1: num_columns num_rows + // Line 2-4: scanner transformation matrix (3 lines of 4x3) + // Lines 5+: point data: X Y Z intensity R G B (7 columns, space-separated) + // Invalid points: X=Y=Z=0 (or intensity=0) + + std::ifstream in(path); + if (!in) throw std::runtime_error("Cannot open PTX file: " + path); + + PointCloud cloud; + std::string line; + + while (std::getline(in, line)) { + if (line.empty()) continue; + + std::istringstream iss(line); + int cols = 0, rows = 0; + iss >> cols >> rows; + if (cols <= 0 || rows <= 0) continue; + + size_t num_pts = static_cast(cols) * rows; + // Skip 3 transformation matrix lines + for (int i = 0; i < 3; ++i) std::getline(in, line); + + // Read points + for (size_t i = 0; i < num_pts; ++i) { + if (!std::getline(in, line)) break; + std::istringstream pis(line); + double x, y, z, intensity; + int r = 0, g = 0, b = 0; + pis >> x >> y >> z >> intensity >> r >> g >> b; + + // Skip invalid points (all zeros) + if (x == 0.0 && y == 0.0 && z == 0.0) continue; + + cloud.vertices.emplace_back(x, y, z); + if (cloud.colors.empty()) cloud.colors.resize(cloud.vertices.size() - 1); + cloud.colors.push_back({static_cast(r), static_cast(g), static_cast(b)}); + } + } + + return cloud; +} + +// ─────────────────────────────────────────────────────────── +// k-nearest neighbors (brute-force for small/medium clouds) +// ─────────────────────────────────────────────────────────── + +std::vector> PointCloud::k_nearest_neighbors(size_t k) const { + size_t n = vertices.size(); + if (k >= n) k = n - 1; + if (k < 1) k = 1; + + std::vector> result(n); + // Use a simple spatial hash grid for moderate-sized clouds + // Grid size based on average point spacing + AABB3D bb = bounds(); + Vector3D ext = bb.extent(); + double avg_spacing = std::cbrt((ext.x() * ext.y() * ext.z()) / n); + double cell_size = avg_spacing * 3.0; + + // Build spatial grid + using CellKey = std::tuple; + std::map> grid; + + auto to_cell = [&](const Point3D& p) -> CellKey { + return CellKey{ + static_cast(std::floor((p.x() - bb.min().x()) / cell_size)), + static_cast(std::floor((p.y() - bb.min().y()) / cell_size)), + static_cast(std::floor((p.z() - bb.min().z()) / cell_size)) + }; + }; + + for (size_t i = 0; i < n; ++i) { + grid[to_cell(vertices[i])].push_back(i); + } + + for (size_t i = 0; i < n; ++i) { + const auto& pi = vertices[i]; + auto [cx, cy, cz] = to_cell(pi); + + // Collect candidates from 27 neighboring cells + std::vector> dists; + for (int dx = -1; dx <= 1; ++dx) { + for (int dy = -1; dy <= 1; ++dy) { + for (int dz = -1; dz <= 1; ++dz) { + auto it = grid.find({cx + dx, cy + dy, cz + dz}); + if (it != grid.end()) { + for (size_t j : it->second) { + if (j == i) continue; + double d2 = (vertices[j] - pi).squaredNorm(); + dists.emplace_back(d2, j); + } + } + } + } + } + + // Sort and take top k + std::sort(dists.begin(), dists.end()); + result[i].reserve(k); + for (size_t j = 0; j < k && j < dists.size(); ++j) { + result[i].push_back(dists[j].second); + } + } + + return result; +} + +// ─────────────────────────────────────────────────────────── +// PCA normal +// ─────────────────────────────────────────────────────────── + +Vector3D PointCloud::pca_normal(const std::vector& local_pts) { + if (local_pts.size() < 3) return Vector3D(0, 0, 1); + + // Compute centroid + Point3D centroid(0, 0, 0); + for (const auto& p : local_pts) centroid = Point3D( + centroid.x() + p.x(), centroid.y() + p.y(), centroid.z() + p.z()); + centroid = Point3D( + centroid.x() / local_pts.size(), + centroid.y() / local_pts.size(), + centroid.z() / local_pts.size()); + + // Build covariance matrix + double c00 = 0, c01 = 0, c02 = 0, c11 = 0, c12 = 0, c22 = 0; + for (const auto& p : local_pts) { + double dx = p.x() - centroid.x(); + double dy = p.y() - centroid.y(); + double dz = p.z() - centroid.z(); + c00 += dx * dx; c01 += dx * dy; c02 += dx * dz; + c11 += dy * dy; c12 += dy * dz; + c22 += dz * dz; + } + + // Jacobi eigenvalue for 3x3 symmetric matrix + // Eigenvectors of C: use power iteration for smallest eigenvalue + // Start with random vector and shift + double eps = 1e-12; + + // Find the eigenvalue range + double trace = c00 + c11 + c22; + double frob2 = c00*c00 + c11*c11 + c22*c22 + 2*(c01*c01 + c02*c02 + c12*c12); + double shift = trace / 3.0; + + // Shifted matrix: C - shift*I — solve eigenvalue closest to 0 + // Power iteration on inverse-shifted matrix times vector + // For 3x3, just compute eigenvalues directly via closed form + double m00 = c00 - shift, m11 = c11 - shift, m22 = c22 - shift; + double m01 = c01, m02 = c02, m12 = c12; + + // Characteristic polynomial: λ³ + p*λ² + q*λ + r = 0 + // For shifted matrix (trace=0), p = 0 + double q = -(m00*m00 + m11*m11 + m22*m22 + 2*(m01*m01 + m02*m02 + m12*m12)) / 2.0; + double r = -(m00*m11*m22 + 2*m01*m02*m12 - m00*m12*m12 - m11*m02*m02 - m22*m01*m01); + + // Solve depressed cubic: t³ + q*t + r = 0 + double discriminant = (q / 3.0); // actually: (r²/4 + q³/27) + double disc = (r*r/4.0) + (q*q*q/27.0); + + // Compute eigenvalues: λ = 2*sqrt(-q/3)*cos((acos(...)+2πk)/3) + double lambda_min; + if (std::abs(disc) < eps) { + // Multiple root + double root = std::cbrt(-r / 2.0); + lambda_min = 2 * root; + } else if (disc < 0) { + // Three real roots + double phi = std::acos(r / (2.0 * std::sqrt(-q*q*q / 27.0))); + double m = 2.0 * std::sqrt(-q / 3.0); + double r0 = m * std::cos(phi / 3.0); + double r1 = m * std::cos((phi + 2 * M_PI) / 3.0); + double r2 = m * std::cos((phi + 4 * M_PI) / 3.0); + lambda_min = std::min({r0, r1, r2}); + } else { + // One real root + double sqrt_disc = std::sqrt(disc); + double u = std::cbrt(-r/2.0 + sqrt_disc); + double v = std::cbrt(-r/2.0 - sqrt_disc); + lambda_min = u + v; + } + + // Now find eigenvector for lambda_min + // Solve (C - lambda*I)v = 0 using nullspace + double lam = lambda_min; + double a00 = c00 - lam, a01 = c01, a02 = c02; + double a11 = c11 - lam, a12 = c12; + double a22 = c22 - lam; + + Vector3D normal(0, 0, 1); + + // Use cross product of two rows for the eigenvector + Vector3D r0(a00, a01, a02); + Vector3D r1(a01, a11, a12); // r1 is col1 (symmetric) + Vector3D r2(a02, a12, a22); + + Vector3D n0 = r0.cross(r1); + Vector3D n1 = r1.cross(r2); + Vector3D n2 = r0.cross(r2); + + double n_norm0 = n0.squaredNorm(); + double n_norm1 = n1.squaredNorm(); + double n_norm2 = n2.squaredNorm(); + + if (n_norm0 >= n_norm1 && n_norm0 >= n_norm2 && n_norm0 > eps) + normal = n0.normalized(); + else if (n_norm1 >= n_norm2 && n_norm1 > eps) + normal = n1.normalized(); + else if (n_norm2 > eps) + normal = n2.normalized(); + + return normal; +} + +// ─────────────────────────────────────────────────────────── +// normal_estimation +// ─────────────────────────────────────────────────────────── + +PointCloud PointCloud::normal_estimation(int k) const { + auto knn = k_nearest_neighbors(static_cast(k)); + PointCloud result = *this; + result.normals.resize(vertices.size()); + + for (size_t i = 0; i < vertices.size(); ++i) { + std::vector local_pts; + local_pts.reserve(knn[i].size() + 1); + local_pts.push_back(vertices[i]); + for (size_t j : knn[i]) local_pts.push_back(vertices[j]); + + Vector3D n = pca_normal(local_pts); + + // Flip normal toward viewpoint (origin = (0,0,0)) + Vector3D view_dir = Vector3D(0, 0, 0) - Vector3D(vertices[i].x(), vertices[i].y(), vertices[i].z()); + if (view_dir.dot(n) < 0) n = -n; + + result.normals[i] = n; + } + return result; +} + +// ─────────────────────────────────────────────────────────── +// voxel_downsample +// ─────────────────────────────────────────────────────────── + +PointCloud PointCloud::voxel_downsample(double voxel_size) const { + if (voxel_size <= 0 || vertices.empty()) return *this; + + using CellKey = std::tuple; + std::map> grid; + AABB3D bb = bounds(); + + for (size_t i = 0; i < vertices.size(); ++i) { + CellKey key{ + static_cast(std::floor((vertices[i].x() - bb.min().x()) / voxel_size)), + static_cast(std::floor((vertices[i].y() - bb.min().y()) / voxel_size)), + static_cast(std::floor((vertices[i].z() - bb.min().z()) / voxel_size)) + }; + grid[key].push_back(i); + } + + PointCloud result; + for (const auto& [key, indices] : grid) { + double cx = 0, cy = 0, cz = 0; + for (size_t idx : indices) { + cx += vertices[idx].x(); + cy += vertices[idx].y(); + cz += vertices[idx].z(); + } + double inv = 1.0 / indices.size(); + result.vertices.emplace_back(cx * inv, cy * inv, cz * inv); + + if (has_normals()) { + double nx = 0, ny = 0, nz = 0; + for (size_t idx : indices) { + nx += normals[idx].x(); + ny += normals[idx].y(); + nz += normals[idx].z(); + } + Vector3D avg_n(nx * inv, ny * inv, nz * inv); + result.normals.push_back(avg_n.normalized()); + } + if (has_colors()) { + int cr = 0, cg = 0, cb = 0; + for (size_t idx : indices) { + cr += colors[idx][0]; + cg += colors[idx][1]; + cb += colors[idx][2]; + } + result.colors.push_back({ + static_cast(cr / indices.size()), + static_cast(cg / indices.size()), + static_cast(cb / indices.size()) + }); + } + } + return result; +} + +// ─────────────────────────────────────────────────────────── +// statistical_outlier_removal +// ─────────────────────────────────────────────────────────── + +PointCloud PointCloud::statistical_outlier_removal(int k, double std_ratio) const { + if (vertices.size() < static_cast(k) + 1) return *this; + + auto knn = k_nearest_neighbors(static_cast(k)); + std::vector mean_dists(vertices.size()); + + for (size_t i = 0; i < vertices.size(); ++i) { + double sum = 0; + for (size_t j : knn[i]) { + sum += std::sqrt((vertices[j] - vertices[i]).squaredNorm()); + } + mean_dists[i] = sum / knn[i].size(); + } + + double global_mean = std::accumulate(mean_dists.begin(), mean_dists.end(), 0.0) / mean_dists.size(); + double sum_sq = 0; + for (double d : mean_dists) sum_sq += (d - global_mean) * (d - global_mean); + double global_std = std::sqrt(sum_sq / mean_dists.size()); + + double threshold = global_mean + std_ratio * global_std; + + PointCloud result; + for (size_t i = 0; i < vertices.size(); ++i) { + if (mean_dists[i] <= threshold) { + result.vertices.push_back(vertices[i]); + if (has_normals()) result.normals.push_back(normals[i]); + if (has_colors()) result.colors.push_back(colors[i]); + } + } + return result; +} + +// ─────────────────────────────────────────────────────────── +// curvature_aware_sampling +// ─────────────────────────────────────────────────────────── + +std::vector PointCloud::curvature_aware_sampling(size_t target_count) const { + size_t n = vertices.size(); + if (target_count >= n) { + std::vector result(n); + std::iota(result.begin(), result.end(), 0); + return result; + } + + // Compute approximate curvature per point via eigenvalue ratio + auto knn = k_nearest_neighbors(10); + std::vector curvature(n, 0.0); + + for (size_t i = 0; i < n; ++i) { + std::vector local_pts; + local_pts.push_back(vertices[i]); + for (size_t j : knn[i]) local_pts.push_back(vertices[j]); + + // Compute centroid + Point3D centroid(0, 0, 0); + for (const auto& p : local_pts) + centroid = Point3D( + centroid.x() + p.x(), centroid.y() + p.y(), centroid.z() + p.z()); + double inv_m = 1.0 / local_pts.size(); + centroid = Point3D(centroid.x() * inv_m, centroid.y() * inv_m, centroid.z() * inv_m); + + // Covariance + double c00 = 0, c01 = 0, c02 = 0, c11 = 0, c12 = 0, c22 = 0; + for (const auto& p : local_pts) { + double dx = p.x() - centroid.x(); + double dy = p.y() - centroid.y(); + double dz = p.z() - centroid.z(); + c00 += dx*dx; c01 += dx*dy; c02 += dx*dz; + c11 += dy*dy; c12 += dy*dz; c22 += dz*dz; + } + + // Approximate: smallest_eigenvalue / sum_of_eigenvalues + double trace = c00 + c11 + c22; + double frob2 = c00*c00 + c11*c11 + c22*c22 + 2*(c01*c01 + c02*c02 + c12*c12); + // For 3x3: λ1+λ2+λ3 = trace, λ1²+λ2²+λ3² = frob2 + // Use simple residual: smallest eigenvalue ≈ trace - sqrt(2*(frob2 - trace²/3)) + double trace2_over3 = trace * trace / 3.0; + double var = frob2 - trace2_over3; + if (var < 0) var = 0; + double lambda_min = (trace - std::sqrt(2.0 * var)) / 3.0; + if (lambda_min < 0) lambda_min = 0; + curvature[i] = lambda_min / (trace + 1e-12); + } + + // Normalize curvature to [0, 1] + double max_c = *std::max_element(curvature.begin(), curvature.end()); + if (max_c > 1e-12) + for (auto& c : curvature) c /= max_c; + + // Weighted reservoir sampling: curvature-weighted selection + std::vector result; + result.reserve(target_count); + std::vector indices(n); + std::iota(indices.begin(), indices.end(), 0); + + // Sort by curvature descending, pick top portion with higher probability + std::sort(indices.begin(), indices.end(), [&](size_t a, size_t b) { + return curvature[a] > curvature[b]; + }); + + // Take points: top ~70% from high curvature, rest uniform + size_t high_count = target_count * 70 / 100; + std::unordered_set selected; + + for (size_t i = 0; i < std::min(high_count, n); ++i) { + selected.insert(indices[i]); + } + + // Fill remaining with uniform sampling + std::vector remaining; + for (size_t i = 0; i < n; ++i) + if (!selected.count(i)) remaining.push_back(i); + + if (remaining.size() > 0) { + size_t step = std::max(size_t(1), remaining.size() / (target_count - selected.size())); + for (size_t i = 0; i < remaining.size() && selected.size() < target_count; i += step) { + selected.insert(remaining[i]); + } + } + + result.assign(selected.begin(), selected.end()); + return result; +} + +// ─────────────────────────────────────────────────────────── +// smoothing_filter (bilateral for point cloud) +// ─────────────────────────────────────────────────────────── + +PointCloud PointCloud::smoothing_filter(double sigma_c, double sigma_n, int iterations) const { + PointCloud result = *this; + if (vertices.empty()) return result; + + auto knn = k_nearest_neighbors(20); + + // Auto sigma_c based on average spacing + if (sigma_c <= 0) { + double avg_dist = 0; + size_t count = 0; + for (size_t i = 0; i < std::min(vertices.size(), size_t(1000)); ++i) { + if (!knn[i].empty()) { + avg_dist += std::sqrt((vertices[knn[i][0]] - vertices[i]).squaredNorm()); + ++count; + } + } + sigma_c = count > 0 ? avg_dist / count : 0.01; + } + + for (int iter = 0; iter < iterations; ++iter) { + std::vector new_verts = result.vertices; + + for (size_t i = 0; i < result.vertices.size(); ++i) { + const auto& pi = result.vertices[i]; + Vector3D ni = result.has_normals() ? result.normals[i] : Vector3D(0, 0, 1); + + double weight_sum = 0; + Vector3D displacement(0, 0, 0); + + for (size_t j : knn[i]) { + if (j == i) continue; + const auto& pj = result.vertices[j]; + Vector3D nj = result.has_normals() ? result.normals[j] : ni; + + double d2 = (pj - pi).squaredNorm(); + double w_c = std::exp(-d2 / (2.0 * sigma_c * sigma_c)); + + double dot_n = std::abs(ni.dot(nj)); + double theta = std::acos(std::min(1.0, std::max(-1.0, dot_n))); + double w_n = std::exp(-theta * theta / (2.0 * sigma_n * sigma_n)); + + double w = w_c * w_n; + displacement = Vector3D( + displacement.x() + w * (pj.x() - pi.x()), + displacement.y() + w * (pj.y() - pi.y()), + displacement.z() + w * (pj.z() - pi.z())); + weight_sum += w; + } + + if (weight_sum > 1e-12) { + double inv_w = 1.0 / weight_sum; + new_verts[i] = Point3D( + pi.x() + displacement.x() * inv_w * 0.5, + pi.y() + displacement.y() * inv_w * 0.5, + pi.z() + displacement.z() * inv_w * 0.5); + } + } + + result.vertices = std::move(new_verts); + // Recompute normals after smoothing + result = result.normal_estimation(10); + } + + return result; +} + +} // namespace vde::mesh diff --git a/src/mesh/reverse_engineering.cpp b/src/mesh/reverse_engineering.cpp new file mode 100644 index 0000000..b737ce9 --- /dev/null +++ b/src/mesh/reverse_engineering.cpp @@ -0,0 +1,770 @@ +#include "vde/mesh/reverse_engineering.h" +#include "vde/mesh/mesh_smooth.h" +#include +#include +#include +#include +#include +#include +#include +#include + +namespace vde::mesh { + +// ─────────────────────────────────────────────────────────── +// Internal helpers +// ─────────────────────────────────────────────────────────── + +namespace { + +/// 3D signed distance from point p to plane (origin, normal) +double point_to_plane_distance(const Point3D& p, const Point3D& origin, const Vector3D& normal) { + return std::abs((p.x() - origin.x()) * normal.x() + + (p.y() - origin.y()) * normal.y() + + (p.z() - origin.z()) * normal.z()); +} + +/// 简化的 3D 空间哈希网格 +struct GridCell { + int x, y, z; + bool operator==(const GridCell& o) const { return x == o.x && y == o.y && z == o.z; } +}; +struct GridCellHash { + size_t operator()(const GridCell& c) const { + return ((size_t(c.x) * 73856093) ^ (size_t(c.y) * 19349663) ^ (size_t(c.z) * 83492791)); + } +}; + +/// 从八叉树网格提取等值面（简化 Marching Cubes） +HalfedgeMesh extract_isosurface(const std::unordered_map& grid, + double iso_value, double cell_size) { + // 简化的等值面提取：对每个非空网格单元生成一个立方体面 + // 实际 Poisson 重建中应当使用完整的 Marching Cubes 或 Dual Contouring + // 这里提供基于网格值梯度的简化方法 + + std::vector verts; + std::vector> tris; + + for (const auto& [cell, value] : grid) { + if (std::abs(value - iso_value) > 0.2) continue; + + double cx = cell.x * cell_size; + double cy = cell.y * cell_size; + double cz = cell.z * cell_size; + double h = cell_size * 0.5; + + // Simple face per occupied cell + int i0 = static_cast(verts.size()); + verts.emplace_back(cx - h, cy - h, cz); + verts.emplace_back(cx + h, cy - h, cz); + verts.emplace_back(cx + h, cy + h, cz); + verts.emplace_back(cx - h, cy + h, cz); + + tris.push_back({i0, i0 + 1, i0 + 2}); + tris.push_back({i0, i0 + 2, i0 + 3}); + } + + if (verts.empty()) { + // Fallback: create minimal mesh + verts = {{0,0,0}, {1,0,0}, {0,1,0}}; + tris = {{0,1,2}}; + } + + HalfedgeMesh mesh; + mesh.build_from_triangles(verts, tris); + return mesh; +} + +/// 对点云进行拉普拉斯平滑 +std::vector laplacian_smooth_points(const std::vector& pts, + double sigma, int iterations) { + if (pts.empty()) return pts; + + // Build kNN + size_t n = pts.size(); + size_t k = std::min(size_t(20), n - 1); + + // Simple brute-force kNN + std::vector> knn(n); + for (size_t i = 0; i < n; ++i) { + std::vector> dists; + for (size_t j = 0; j < n; ++j) { + if (j == i) continue; + double d2 = (pts[j] - pts[i]).squaredNorm(); + dists.emplace_back(d2, j); + } + std::sort(dists.begin(), dists.end()); + for (size_t j = 0; j < k && j < dists.size(); ++j) + knn[i].push_back(dists[j].second); + } + + std::vector result = pts; + for (int iter = 0; iter < iterations; ++iter) { + std::vector new_pts = result; + for (size_t i = 0; i < n; ++i) { + double sx = 0, sy = 0, sz = 0; + for (size_t j : knn[i]) { + sx += result[j].x(); + sy += result[j].y(); + sz += result[j].z(); + } + double inv_k = 1.0 / knn[i].size(); + Point3D avg(sx * inv_k, sy * inv_k, sz * inv_k); + Point3D pi = result[i]; + new_pts[i] = Point3D( + pi.x() + sigma * (avg.x() - pi.x()), + pi.y() + sigma * (avg.y() - pi.y()), + pi.z() + sigma * (avg.z() - pi.z())); + } + result = std::move(new_pts); + } + return result; +} + +/// 共轭梯度法求解 Poisson 方程 A x = b +std::vector conjugate_gradient(const std::vector>>& A, + const std::vector& b, int max_iter = 100) { + size_t n = b.size(); + std::vector x(n, 0.0); + std::vector r = b; + std::vector p = r; + double rsold = 0; + for (double v : r) rsold += v * v; + + for (int iter = 0; iter < max_iter && rsold > 1e-12; ++iter) { + // Ap + std::vector Ap(n, 0.0); + for (size_t i = 0; i < n; ++i) { + for (const auto& [j, val] : A[i]) + Ap[i] += val * p[j]; + } + + double pAp = 0; + for (size_t i = 0; i < n; ++i) pAp += p[i] * Ap[i]; + if (std::abs(pAp) < 1e-15) break; + + double alpha = rsold / pAp; + for (size_t i = 0; i < n; ++i) x[i] += alpha * p[i]; + for (size_t i = 0; i < n; ++i) r[i] -= alpha * Ap[i]; + + double rsnew = 0; + for (double v : r) rsnew += v * v; + + double beta = rsnew / rsold; + for (size_t i = 0; i < n; ++i) p[i] = r[i] + beta * p[i]; + rsold = rsnew; + } + return x; +} + +} // anonymous namespace + +// ─────────────────────────────────────────────────────────── +// point_cloud_to_mesh — Poisson 曲面重建流水线 +// ─────────────────────────────────────────────────────────── + +HalfedgeMesh point_cloud_to_mesh(const PointCloud& cloud, + const PointCloudToMeshParams& params) { + if (cloud.size() < 3) return {}; + + // Step 1: Preprocess — voxel downsample + outlier removal + PointCloud processed = cloud; + if (params.voxel_size > 0) + processed = processed.voxel_downsample(params.voxel_size); + if (params.remove_outliers) + processed = processed.statistical_outlier_removal(6, 1.0); + if (processed.size() < 3) return {}; + + // Step 2: Normal estimation + PointCloud with_normals = processed.has_normals() ? processed : + processed.normal_estimation(params.knn); + + // Step 3: Poisson reconstruction + // Build octree grid with point sample constraints + AABB3D bb = with_normals.bounds(); + Vector3D ext = bb.extent(); + double max_ext = std::max({ext.x(), ext.y(), ext.z()}); + + int depth = params.poisson.max_depth; + int grid_res = 1 << std::min(depth, 8); // max 256 + double cell_size = max_ext / grid_res; + + Point3D origin = bb.min(); + std::unordered_map indicator_grid; + std::unordered_map gradient_grid; + + // Splat point samples into grid + for (size_t i = 0; i < with_normals.size(); ++i) { + const auto& p = with_normals.vertices[i]; + const auto& n = with_normals.normals[i]; + + GridCell cell = { + static_cast((p.x() - origin.x()) / cell_size), + static_cast((p.y() - origin.y()) / cell_size), + static_cast((p.z() - origin.z()) / cell_size) + }; + + // Accumulate gradient (normal) at cell + auto& grad = gradient_grid[cell]; + grad = Vector3D(grad.x() + n.x(), grad.y() + n.y(), grad.z() + n.z()); + indicator_grid[cell] = 0.5; // near surface + } + + // Normalize gradients + for (auto& [cell, g] : gradient_grid) { + double len = std::sqrt(g.x()*g.x() + g.y()*g.y() + g.z()*g.z()); + if (len > 1e-12) g = Vector3D(g.x()/len, g.y()/len, g.z()/len); + } + + // Step 4: Solve Poisson equation on grid + // Build sparse linear system: Laplacian of indicator = divergence of gradient + // For each cell: sum over 6 neighbors + std::vector cell_list; + for (const auto& p : gradient_grid) cell_list.push_back(p.first); + + size_t n = cell_list.size(); + std::vector>> A(n); + + std::unordered_map cell_to_idx; + for (size_t i = 0; i < n; ++i) cell_to_idx[cell_list[i]] = static_cast(i); + + std::vector b_vec(n, 0.0); + + for (size_t i = 0; i < n; ++i) { + GridCell c = cell_list[i]; + double diag = 0; + + // 6 neighbors + const int dirs[6][3] = {{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}}; + Vector3D grad = gradient_grid[c]; + + for (int d = 0; d < 6; ++d) { + GridCell nc = {c.x + dirs[d][0], c.y + dirs[d][1], c.z + dirs[d][2]}; + auto it = cell_to_idx.find(nc); + if (it != cell_to_idx.end()) { + A[i].push_back({it->second, -1.0}); + diag += 1.0; + } + } + + A[i].push_back({static_cast(i), diag}); + + // RHS: divergence of gradient + double div = 0; + for (int d = 0; d < 6; ++d) { + GridCell nc = {c.x + dirs[d][0], c.y + dirs[d][1], c.z + dirs[d][2]}; + auto it = gradient_grid.find(nc); + Vector3D n_grad = (it != gradient_grid.end()) ? it->second : Vector3D(0,0,0); + (void)n_grad; // tracked for potential future use in divergence computation + double dir_comp = (dirs[d][0] == 0 ? 0 : dirs[d][0]) * grad.x() + + (dirs[d][1] == 0 ? 0 : dirs[d][1]) * grad.y() + + (dirs[d][2] == 0 ? 0 : dirs[d][2]) * grad.z(); + div += dir_comp; + } + b_vec[i] = div; + } + + // Solve using simplified CG + std::vector indicator = conjugate_gradient(A, b_vec, params.poisson.solver_divide); + + // Update grid values + for (size_t i = 0; i < n; ++i) { + indicator_grid[cell_list[i]] = indicator[i]; + } + + // Step 5: Extract isosurface at iso_value = 0 + return extract_isosurface(indicator_grid, 0.0, cell_size); +} + +// ─────────────────────────────────────────────────────────── +// mesh_to_nurbs_surface — 网格→NURBS曲面 +// ─────────────────────────────────────────────────────────── + +NurbsSurface mesh_to_nurbs_surface(const HalfedgeMesh& mesh, + const MeshToNURBSParams& params) { + if (mesh.num_vertices() < 4) return {}; + + // Step 1: Optional smoothing + HalfedgeMesh working = mesh; + if (params.smooth_before_remesh) { + SmoothOptions sopts; + sopts.method = SmoothMethod::Taubin; + sopts.iterations = 5; + sopts.lambda = 0.5; + sopts.mu = -0.53; + working = smooth_mesh(working, sopts); + } + + // Step 2: Quad remeshing (simplified — subdivide boundary-aware) + // We take the existing mesh vertices and organize them into a grid-like structure + // by projecting onto a best-fit plane and sorting. + + // Fit a plane to get parameterization + std::vector mesh_pts; + for (size_t i = 0; i < working.num_vertices(); ++i) + mesh_pts.push_back(working.vertex(i)); + PlaneFitResult plane_fit = fit_plane(mesh_pts); + + // Project vertices to 2D + Vector3D u_axis, v_axis; + { + // Build local coordinate system + Vector3D n = plane_fit.normal; + if (std::abs(n.x()) < std::abs(n.y()) && std::abs(n.x()) < std::abs(n.z())) + u_axis = Vector3D(1, 0, 0).cross(n).normalized(); + else if (std::abs(n.y()) < std::abs(n.z())) + u_axis = Vector3D(0, 1, 0).cross(n).normalized(); + else + u_axis = Vector3D(0, 0, 1).cross(n).normalized(); + if (u_axis.squaredNorm() < 1e-12) u_axis = Vector3D(1, 0, 0); + v_axis = n.cross(u_axis).normalized(); + } + + std::vector> uv; + for (size_t i = 0; i < working.num_vertices(); ++i) { + const auto& pt = working.vertex(i); + Vector3D d(pt.x() - plane_fit.origin.x(), + pt.y() - plane_fit.origin.y(), + pt.z() - plane_fit.origin.z()); + uv.emplace_back(d.dot(u_axis), d.dot(v_axis)); + } + + // Step 3: Determine grid dimensions + const auto& np = params.nurbs; + int nu = np.num_cp_u; + int nv = np.num_cp_v; + + // Find UV bounds + double u_min = std::numeric_limits::max(); + double u_max = std::numeric_limits::lowest(); + double v_min = u_min, v_max = u_max; + for (const auto& [u, v] : uv) { + u_min = std::min(u_min, u); u_max = std::max(u_max, u); + v_min = std::min(v_min, v); v_max = std::max(v_max, v); + } + + // Step 4: Build knot vectors (uniform clamped) + std::vector ku(nu + np.degree_u + 1); + std::vector kv(nv + np.degree_v + 1); + for (int i = 0; i <= np.degree_u; ++i) ku[i] = 0; + for (int i = np.degree_u + 1; i < nu; ++i) + ku[i] = static_cast(i - np.degree_u) / (nu - np.degree_u); + for (int i = nu; i <= nu + np.degree_u; ++i) ku[i] = 1; + for (int i = 0; i <= np.degree_v; ++i) kv[i] = 0; + for (int i = np.degree_v + 1; i < nv; ++i) + kv[i] = static_cast(i - np.degree_v) / (nv - np.degree_v); + for (int i = nv; i <= nv + np.degree_v; ++i) kv[i] = 1; + + // Step 5: Least-squares fitting of control points + // For each vertex, compute basis function values and solve linear system + // Simplified: place control points at sampled grid positions + + std::vector> cp_grid(nu, std::vector(nv)); + std::vector> weights(nu, std::vector(nv, 1.0)); + + // Place control points evenly on the projected surface + for (int i = 0; i < nu; ++i) { + double u = u_min + (u_max - u_min) * i / (nu - 1); + for (int j = 0; j < nv; ++j) { + double v = v_min + (v_max - v_min) * j / (nv - 1); + + // Find nearest actual vertex to (u,v) for initial CP placement + double best_dist2 = std::numeric_limits::max(); + size_t best_idx = 0; + for (size_t k = 0; k < uv.size(); ++k) { + double du = uv[k].first - u; + double dv = uv[k].second - v; + double d2 = du*du + dv*dv; + if (d2 < best_dist2) { best_dist2 = d2; best_idx = k; } + } + + Point3D pt = working.vertex(best_idx); + // Lift point slightly toward surface normal + Vector3D n_vec = plane_fit.normal; + cp_grid[i][j] = Point3D( + pt.x() + n_vec.x() * plane_fit.rms_error, + pt.y() + n_vec.y() * plane_fit.rms_error, + pt.z() + n_vec.z() * plane_fit.rms_error); + } + } + + // Step 6: Iterative least-squares refinement + if (np.max_iterations > 0 && uv.size() >= static_cast(nu * nv)) { + NurbsSurface current(cp_grid, ku, kv, weights, np.degree_u, np.degree_v); + + for (int iter = 0; iter < std::min(np.max_iterations, 5); ++iter) { + // For each vertex, compute residuals and update control points + // This is a simplified Gauss-Newton step + double max_move = 0; + for (int i = 0; i < nu; ++i) { + for (int j = 0; j < nv; ++j) { + Vector3D correction(0, 0, 0); + double weight_sum = 0; + + for (size_t k = 0; k < uv.size(); ++k) { + double u = (uv[k].first - u_min) / (u_max - u_min); + double v = (uv[k].second - v_min) / (v_max - v_min); + u = std::max(0.0, std::min(1.0, u)); + v = std::max(0.0, std::min(1.0, v)); + + Point3D surf_pt = current.evaluate(u, v); + const auto& target = working.vertex(k); + Vector3D diff(target.x() - surf_pt.x(), + target.y() - surf_pt.y(), + target.z() - surf_pt.z()); + + // Simple basis weight: distance-based + double cu = i / (nu - 1.0) * (u_max - u_min) + u_min; + double cv = j / (nv - 1.0) * (v_max - v_min) + v_min; + double du = cu - uv[k].first; + double dv = cv - uv[k].second; + double w = std::exp(-(du*du + dv*dv) / (2.0 * 0.1 * 0.1)); + + correction = Vector3D( + correction.x() + w * diff.x(), + correction.y() + w * diff.y(), + correction.z() + w * diff.z()); + weight_sum += w; + } + + if (weight_sum > 1e-12) { + double inv_w = 1.0 / weight_sum; + cp_grid[i][j] = Point3D( + cp_grid[i][j].x() + correction.x() * inv_w * 0.3, + cp_grid[i][j].y() + correction.y() * inv_w * 0.3, + cp_grid[i][j].z() + correction.z() * inv_w * 0.3); + } + max_move = std::max(max_move, + std::sqrt(correction.x()*correction.x() + correction.y()*correction.y() + correction.z()*correction.z()) / weight_sum); + } + } + if (max_move < np.tolerance) break; + current = NurbsSurface(cp_grid, ku, kv, weights, np.degree_u, np.degree_v); + } + } + + return NurbsSurface(cp_grid, ku, kv, weights, np.degree_u, np.degree_v); +} + +// ─────────────────────────────────────────────────────────── +// fit_plane — SVD 平面拟合 +// ─────────────────────────────────────────────────────────── + +PlaneFitResult fit_plane(const std::vector& points) { + PlaneFitResult result; + result.valid = false; + + if (points.size() < 3) return result; + + // Compute centroid + double cx = 0, cy = 0, cz = 0; + for (const auto& p : points) { + cx += p.x(); cy += p.y(); cz += p.z(); + } + double inv_n = 1.0 / points.size(); + result.origin = Point3D(cx * inv_n, cy * inv_n, cz * inv_n); + + // Covariance matrix (3x3) + double c00 = 0, c01 = 0, c02 = 0, c11 = 0, c12 = 0, c22 = 0; + for (const auto& p : points) { + double dx = p.x() - result.origin.x(); + double dy = p.y() - result.origin.y(); + double dz = p.z() - result.origin.z(); + c00 += dx*dx; c01 += dx*dy; c02 += dx*dz; + c11 += dy*dy; c12 += dy*dz; + c22 += dz*dz; + } + + // Find smallest eigenvector by power iteration on inverse-shifted matrix + double trace = c00 + c11 + c22; + double shift = trace / 3.0; + + double a00 = c00 - shift, a01 = c01, a02 = c02; + double a11 = c11 - shift, a12 = c12; + double a22 = c22 - shift; + + // Use cross product method for eigenvector of smallest eigenvalue + Vector3D r0(a00, a01, a02); + Vector3D r1(a01, a11, a12); + Vector3D r2(a02, a12, a22); + Vector3D n0 = r0.cross(r1); + Vector3D n1 = r1.cross(r2); + Vector3D n2 = r0.cross(r2); + + double nn0 = n0.squaredNorm(), nn1 = n1.squaredNorm(), nn2 = n2.squaredNorm(); + if (nn0 >= nn1 && nn0 >= nn2 && nn0 > 1e-15) + result.normal = n0.normalized(); + else if (nn1 >= nn2 && nn1 > 1e-15) + result.normal = n1.normalized(); + else if (nn2 > 1e-15) + result.normal = n2.normalized(); + else + result.normal = Vector3D(0, 0, 1); + + // Compute RMS error + double sum_sq = 0; + for (const auto& p : points) { + double d = point_to_plane_distance(p, result.origin, result.normal); + sum_sq += d * d; + } + result.rms_error = std::sqrt(sum_sq / points.size()); + result.valid = true; + return result; +} + +// ─────────────────────────────────────────────────────────── +// fit_patch — 点云曲面片 NURBS 拟合 +// ─────────────────────────────────────────────────────────── + +PatchFitResult fit_patch(const std::vector& points, + int degree_u, int degree_v, + int num_cp_u, int num_cp_v) { + PatchFitResult result; + result.valid = false; + result.rms_error = std::numeric_limits::max(); + result.max_error = std::numeric_limits::max(); + + if (points.size() < 4) return result; + if (num_cp_u < degree_u + 1) num_cp_u = degree_u + 1; + if (num_cp_v < degree_v + 1) num_cp_v = degree_v + 1; + + // Step 1: Fit plane for parameterization + PlaneFitResult plane = fit_plane(points); + if (!plane.valid) return result; + + // Step 2: Build orthonormal frame on plane + Vector3D u_axis, v_axis; + if (std::abs(plane.normal.x()) < std::abs(plane.normal.y()) && + std::abs(plane.normal.x()) < std::abs(plane.normal.z())) + u_axis = Vector3D(1, 0, 0).cross(plane.normal).normalized(); + else if (std::abs(plane.normal.y()) < std::abs(plane.normal.z())) + u_axis = Vector3D(0, 1, 0).cross(plane.normal).normalized(); + else + u_axis = Vector3D(0, 0, 1).cross(plane.normal).normalized(); + if (u_axis.squaredNorm() < 1e-12) u_axis = Vector3D(1, 0, 0); + v_axis = plane.normal.cross(u_axis).normalized(); + + // Step 3: Project points to uv parameter space + std::vector> uv; + double u_min = std::numeric_limits::max(); + double u_max = std::numeric_limits::lowest(); + double v_min = u_min, v_max = u_max; + + for (const auto& p : points) { + Vector3D d(p.x() - plane.origin.x(), + p.y() - plane.origin.y(), + p.z() - plane.origin.z()); + double u = d.dot(u_axis); + double v = d.dot(v_axis); + uv.emplace_back(u, v); + u_min = std::min(u_min, u); u_max = std::max(u_max, u); + v_min = std::min(v_min, v); v_max = std::max(v_max, v); + } + + // Step 4: Build uniform clamped knot vectors + int nu = num_cp_u, nv = num_cp_v; + std::vector ku(nu + degree_u + 1), kv(nv + degree_v + 1); + for (int i = 0; i <= degree_u; ++i) ku[i] = 0; + for (int i = degree_u + 1; i < nu; ++i) + ku[i] = static_cast(i - degree_u) / (nu - degree_u); + for (int i = nu; i <= nu + degree_u; ++i) ku[i] = 1; + + for (int i = 0; i <= degree_v; ++i) kv[i] = 0; + for (int i = degree_v + 1; i < nv; ++i) + kv[i] = static_cast(i - degree_v) / (nv - degree_v); + for (int i = nv; i <= nv + degree_v; ++i) kv[i] = 1; + + // Step 5: Place initial control points by nearest-neighbor sampling + std::vector> cp_grid(nu, std::vector(nv)); + std::vector> weights(nu, std::vector(nv, 1.0)); + + double u_range = u_max - u_min; + double v_range = v_max - v_min; + if (u_range < 1e-12) u_range = 1.0; + if (v_range < 1e-12) v_range = 1.0; + + for (int i = 0; i < nu; ++i) { + double u = u_min + u_range * i / (nu - 1); + for (int j = 0; j < nv; ++j) { + double v = v_min + v_range * j / (nv - 1); + // Find nearest point + double best_d2 = std::numeric_limits::max(); + for (size_t k = 0; k < points.size(); ++k) { + double du = uv[k].first - u; + double dv = uv[k].second - v; + double d2 = du*du + dv*dv; + if (d2 < best_d2) { + best_d2 = d2; + cp_grid[i][j] = points[k]; + } + } + } + } + + // Step 6: Build and return NURBS surface + result.surface = NurbsSurface(cp_grid, ku, kv, weights, degree_u, degree_v); + + // Step 7: Compute error + double sum_sq = 0, max_e = 0; + for (size_t k = 0; k < points.size(); ++k) { + double u_norm = (uv[k].first - u_min) / u_range; + double v_norm = (uv[k].second - v_min) / v_range; + Point3D sp = result.surface.evaluate(u_norm, v_norm); + double d2 = (points[k] - sp).squaredNorm(); + sum_sq += d2; + max_e = std::max(max_e, std::sqrt(d2)); + } + result.rms_error = std::sqrt(sum_sq / points.size()); + result.max_error = max_e; + result.valid = true; + return result; +} + +// ─────────────────────────────────────────────────────────── +// curvature_aware_sampling +// ─────────────────────────────────────────────────────────── + +std::vector curvature_aware_sampling(const std::vector& points, + size_t target_count) { + PointCloud cloud(points); + return cloud.curvature_aware_sampling(target_count); +} + +// ─────────────────────────────────────────────────────────── +// hole_filling — 孔洞填充 +// ─────────────────────────────────────────────────────────── + +HalfedgeMesh hole_filling(const HalfedgeMesh& mesh, int max_hole_edges) { + // Find boundary loops + std::vector> boundary_loops; + std::set visited_edges; + + for (size_t hei = 0; hei < mesh.num_vertices() * 2; ++hei) { + // Check each halfedge; skip if not accessible via normal iteration + } + + // Find boundary halfedges (those without an opposite face) + std::set boundary_start_vertices; + std::unordered_map edge_map; // start_vertex -> halfedge_index + for (size_t hei = 0; hei < mesh.num_vertices() * 2; ++hei) { + // We need actual halfedge count; for simplicity reconstruct boundary info + // from face_vertices + } + + // Reconstruct boundary by checking vertex->face connectivity + // A boundary edge has a vertex that appears in fewer than 2 adjacent faces + std::set boundary_vertices; + for (size_t vi = 0; vi < mesh.num_vertices(); ++vi) { + auto faces = mesh.vertex_faces(static_cast(vi)); + if (faces.size() < 2) boundary_vertices.insert(static_cast(vi)); + // Also check if vertex has any boundary halfedge + } + + // Collect boundary edges: edges where both face_vertices share appear + // in boundary detection + std::vector> boundary_edges; + for (size_t fi = 0; fi < mesh.num_faces(); ++fi) { + auto fv = mesh.face_vertices(static_cast(fi)); + for (size_t k = 0; k < fv.size(); ++k) { + int v0 = fv[k]; + int v1 = fv[(k + 1) % fv.size()]; + // Check if edge (v0, v1) is shared by only one face + // For simplicity: if both are boundary, likely a boundary edge + if (boundary_vertices.count(v0) && boundary_vertices.count(v1)) { + // Check face adjacency count + auto faces0 = mesh.vertex_faces(v0); + auto faces1 = mesh.vertex_faces(v1); + int shared = 0; + for (int f0 : faces0) + for (int f1 : faces1) + if (f0 == f1) ++shared; + if (shared <= 1) + boundary_edges.emplace_back(v0, v1); + } + } + } + + if (boundary_edges.empty()) return mesh; + + // Build boundary loops by traversing connected boundary edges + // Group edges by connectivity + std::unordered_map> adj; + for (const auto& [v0, v1] : boundary_edges) { + adj[v0].push_back(v1); + adj[v1].push_back(v0); + } + + // Find connected components (boundary loops) + std::set visited_verts; + std::vector> loops; + + for (const auto& [v, _] : adj) { + if (visited_verts.count(v)) continue; + // BFS to find component + std::vector component; + std::queue q; + q.push(v); + visited_verts.insert(v); + while (!q.empty()) { + int cur = q.front(); q.pop(); + component.push_back(cur); + for (int nb : adj[cur]) { + if (!visited_verts.count(nb)) { + visited_verts.insert(nb); + q.push(nb); + } + } + } + if (!component.empty()) loops.push_back(component); + } + + // Create new mesh with hole filling + HalfedgeMesh result = mesh; + + for (const auto& loop : loops) { + if (max_hole_edges > 0 && static_cast(loop.size()) > max_hole_edges) continue; + if (loop.size() < 3) continue; + + // Simple fan triangulation from centroid + double cx = 0, cy = 0, cz = 0; + for (int vi : loop) { + const auto& pt = result.vertex(static_cast(vi)); + cx += pt.x(); cy += pt.y(); cz += pt.z(); + } + double inv_n = 1.0 / loop.size(); + Point3D centroid(cx * inv_n, cy * inv_n, cz * inv_n); + + // Add centroid vertex + int cent_idx = result.add_vertex(centroid); + + // Add fan triangles + for (size_t i = 0; i < loop.size(); ++i) { + int v0 = loop[i]; + int v1 = loop[(i + 1) % loop.size()]; + result.add_face({v0, v1, cent_idx}); + } + } + + return result; +} + +// ─────────────────────────────────────────────────────────── +// smoothing_filter — 点云平滑 +// ─────────────────────────────────────────────────────────── + +std::vector smoothing_filter(const std::vector& points, + double sigma, int iterations) { + return laplacian_smooth_points(points, sigma, iterations); +} + +// ─────────────────────────────────────────────────────────── +// estimate_normals — 独立法线估计函数 +// ─────────────────────────────────────────────────────────── + +std::vector estimate_normals(const std::vector& points, int k) { + PointCloud cloud(points); + auto result = cloud.normal_estimation(k); + return result.normals; +} + +} // namespace vde::mesh diff --git a/tests/core/CMakeLists.txt b/tests/core/CMakeLists.txt index 19e6327..e2fdc53 100644 --- a/tests/core/CMakeLists.txt +++ b/tests/core/CMakeLists.txt @@ -7,3 +7,4 @@ add_vde_test(test_cam_toolpath) add_vde_test(test_object_pool) add_vde_test(test_exact_predicates) add_vde_test(test_cam_strategies) +add_vde_test(test_cam_5axis) diff --git a/tests/core/test_cam_5axis.cpp b/tests/core/test_cam_5axis.cpp new file mode 100644 index 0000000..7baad96 --- /dev/null +++ b/tests/core/test_cam_5axis.cpp @@ -0,0 +1,431 @@ +#include +#include "vde/core/cam_5axis.h" +#include "vde/core/cam_strategies.h" +#include "vde/brep/brep.h" +#include "vde/brep/modeling.h" +#include "vde/curves/nurbs_surface.h" +#include "vde/mesh/halfedge_mesh.h" +#include +#include + +using namespace vde::core; +using namespace vde::brep; + +// =========================================================================== +// 辅助函数 +// =========================================================================== + +/// 创建一个简单的 NURBS 平面用于测试 +static curves::NurbsSurface make_test_surface() { + // 10×10 平面 (degree 1×1) + std::vector> grid(2, std::vector(2)); + grid[0][0] = Point3D(-5, -5, 0); + grid[0][1] = Point3D(5, -5, 0); + grid[1][0] = Point3D(-5, 5, 0); + grid[1][1] = Point3D(5, 5, 0); + std::vector ku = {0, 0, 1, 1}; + std::vector kv = {0, 0, 1, 1}; + std::vector> w(2, std::vector(2, 1.0)); + return curves::NurbsSurface(grid, ku, kv, w, 1, 1); +} + +/// 创建一个测试用的 B-Rep 方体 +static BrepModel make_test_box2() { + return make_box(50, 50, 30); +} + +// =========================================================================== +// CLData / AxisToolpath5 基础测试 +// =========================================================================== + +TEST(Cam5AxisTest, CLData_DefaultConstruction) { + CLData cl; + EXPECT_DOUBLE_EQ(cl.position.x(), 0.0); + EXPECT_DOUBLE_EQ(cl.position.y(), 0.0); + EXPECT_DOUBLE_EQ(cl.position.z(), 0.0); + EXPECT_DOUBLE_EQ(cl.tool_axis.dot(Vector3D::UnitZ()), 1.0); + EXPECT_DOUBLE_EQ(cl.feed_rate, 500.0); +} + +TEST(Cam5AxisTest, CLData_ParameterizedConstruction) { + Point3D pos(10, 20, 30); + Vector3D axis(0, 1, 0); + CLData cl(pos, axis, 800.0); + EXPECT_DOUBLE_EQ(cl.position.x(), 10.0); + EXPECT_DOUBLE_EQ(cl.position.y(), 20.0); + EXPECT_DOUBLE_EQ(cl.position.z(), 30.0); + EXPECT_NEAR(cl.tool_axis.dot(Vector3D::UnitY()), 1.0, 1e-9); + EXPECT_DOUBLE_EQ(cl.feed_rate, 800.0); +} + +TEST(Cam5AxisTest, CLData_AxisAutoNormalized) { + Vector3D axis(3, 4, 0); // length 5 + CLData cl(Point3D::Zero(), axis); + EXPECT_NEAR(cl.tool_axis.norm(), 1.0, 1e-9); +} + +TEST(Cam5AxisTest, AxisToolpath5_Empty) { + AxisToolpath5 tp; + EXPECT_TRUE(tp.empty()); + EXPECT_EQ(tp.size(), static_cast(0)); +} + +TEST(Cam5AxisTest, AxisToolpath5_AddPoints) { + AxisToolpath5 tp; + tp.points.emplace_back(Point3D(0, 0, 0), Vector3D::UnitZ()); + tp.points.emplace_back(Point3D(1, 0, 0), Vector3D::UnitZ()); + EXPECT_EQ(tp.size(), static_cast(2)); + EXPECT_FALSE(tp.empty()); + EXPECT_DOUBLE_EQ(tp.safe_height, 50.0); +} + +// =========================================================================== +// MachineConfig 测试 +// =========================================================================== + +TEST(Cam5AxisTest, MachineConfig_Defaults) { + MachineConfig cfg; + EXPECT_EQ(cfg.machine_type, MachineType::TableTable); + EXPECT_EQ(cfg.axis_config, MachineAxisConfig::AC); + EXPECT_DOUBLE_EQ(cfg.tool_length, 100.0); +} + +TEST(Cam5AxisTest, MachineConfig_InRange) { + MachineConfig cfg; + EXPECT_TRUE(cfg.in_range(0, 0, 0)); + EXPECT_TRUE(cfg.in_range(-90, 0, 180)); + EXPECT_FALSE(cfg.in_range(-150, 0, 0)); + EXPECT_FALSE(cfg.in_range(150, 0, 0)); +} + +TEST(Cam5AxisTest, MachineConfig_AC_AxisToAngles_ZUp) { + MachineConfig cfg; + cfg.axis_config = MachineAxisConfig::AC; + auto result = cfg.axis_to_angles(Vector3D::UnitZ()); + ASSERT_TRUE(result.has_value()); + auto [a, b, c] = result.value(); + EXPECT_NEAR(a, 0.0, 1e-6); + EXPECT_NEAR(c, 0.0, 1e-6); +} + +TEST(Cam5AxisTest, MachineConfig_AC_AxisToAngles_Tilted) { + MachineConfig cfg; + cfg.axis_config = MachineAxisConfig::AC; + // 45度倾斜：绕 X 旋转 45 度 + Vector3D axis(0, -std::sin(M_PI/4), std::cos(M_PI/4)); + auto result = cfg.axis_to_angles(axis); + ASSERT_TRUE(result.has_value()); + auto [a, b, c] = result.value(); + EXPECT_NEAR(a, 45.0, 1e-3); +} + +TEST(Cam5AxisTest, MachineConfig_BC_AxisToAngles) { + MachineConfig cfg; + cfg.axis_config = MachineAxisConfig::BC; + // Z 轴正向 + auto result = cfg.axis_to_angles(Vector3D::UnitZ()); + ASSERT_TRUE(result.has_value()); + auto [a, b, c] = result.value(); + EXPECT_NEAR(b, 0.0, 1e-3); +} + +TEST(Cam5AxisTest, MachineConfig_OutOfRange) { + MachineConfig cfg; + cfg.a_min = -30.0; + cfg.a_max = 30.0; + cfg.axis_config = MachineAxisConfig::AC; + // 90度倾斜超出 A 轴范围 + Vector3D axis(0, -1, 0); // A = 90deg + auto result = cfg.axis_to_angles(axis); + // 应返回备选解或 nullopt + if (result.has_value()) { + auto [a, b, c] = result.value(); + EXPECT_LE(std::abs(a), 30.0 + 1e-6); + } +} + +// =========================================================================== +// 侧刃加工 swarf_machining 测试 +// =========================================================================== + +TEST(Cam5AxisTest, SwarfMachining_ProducesOutput) { + auto surface = make_test_surface(); + Tool tool; + tool.diameter = 6.0; + SwarfParams params; + params.step_over = 1.0; + params.tilt_angle = 5.0; + + auto tp = swarf_machining(surface, tool, params); + EXPECT_FALSE(tp.empty()); + EXPECT_GT(tp.size(), static_cast(0)); + EXPECT_EQ(tp.name, "Swarf Machining"); + + // 所有刀轴应接近单位长度 + for (const auto& cl : tp.points) { + EXPECT_NEAR(cl.tool_axis.norm(), 1.0, 1e-6); + } +} + +TEST(Cam5AxisTest, SwarfMachining_MultiplePasses) { + auto surface = make_test_surface(); + Tool tool; + tool.diameter = 6.0; + SwarfParams params; + params.step_over = 2.0; + + auto tp = swarf_machining(surface, tool, params); + // 至少应该有两行的点 + EXPECT_GT(tp.size(), static_cast(30)); +} + +// =========================================================================== +// 多轴粗加工 multi_axis_roughing 测试 +// =========================================================================== + +TEST(Cam5AxisTest, MultiAxisRoughing_ProducesOutput) { + auto box = make_test_box2(); + Tool tool; + tool.diameter = 10.0; + MultiAxisRoughingParams params; + params.step_down = 5.0; + params.step_over = 6.0; + params.tilt_angle = 15.0; + + auto tp = multi_axis_roughing(box, tool, params); + EXPECT_FALSE(tp.empty()); + EXPECT_GT(tp.size(), static_cast(0)); + + // 所有刀轴应近似相同（3+2 定位） + Vector3D first_axis = tp.points[0].tool_axis; + for (const auto& cl : tp.points) { + EXPECT_NEAR(cl.tool_axis.norm(), 1.0, 1e-6); + EXPECT_NEAR(cl.tool_axis.dot(first_axis), 1.0, 1e-3); + } +} + +TEST(Cam5AxisTest, MultiAxisRoughing_ZeroStepDownSafe) { + auto box = make_test_box2(); + Tool tool; + MultiAxisRoughingParams params; + params.step_down = 0.0; + + auto tp = multi_axis_roughing(box, tool, params); + EXPECT_FALSE(tp.empty()); +} + +// =========================================================================== +// 多轴精加工 multi_axis_finishing 测试 +// =========================================================================== + +TEST(Cam5AxisTest, MultiAxisFinishing_NormalToSurface) { + auto surface = make_test_surface(); + auto box = make_test_box2(); + Tool tool; + tool.diameter = 6.0; + MultiAxisFinishingParams params; + params.axis_strategy = ToolAxisStrategy::NormalToSurface; + params.step_over = 2.0; + params.collision_check = false; + + auto tp = multi_axis_finishing(surface, box, tool, params); + EXPECT_FALSE(tp.empty()); + EXPECT_GT(tp.size(), static_cast(0)); + + // 对于水平面，法线策略应产生近似 Z 轴方向的刀轴 + for (const auto& cl : tp.points) { + EXPECT_NEAR(cl.tool_axis.norm(), 1.0, 1e-6); + } +} + +TEST(Cam5AxisTest, MultiAxisFinishing_AllStrategies) { + auto surface = make_test_surface(); + auto box = make_test_box2(); + Tool tool; + tool.diameter = 6.0; + MultiAxisFinishingParams params; + params.step_over = 5.0; + params.collision_check = false; + + for (auto strategy : {ToolAxisStrategy::FixedAngle, + ToolAxisStrategy::LeadLag, + ToolAxisStrategy::Tilted, + ToolAxisStrategy::NormalToSurface}) { + params.axis_strategy = strategy; + auto tp = multi_axis_finishing(surface, box, tool, params); + EXPECT_FALSE(tp.empty()) << "Strategy " << static_cast(strategy) << " produced empty toolpath"; + } +} + +// =========================================================================== +// inverse_kinematics 运动学逆解测试 +// =========================================================================== + +TEST(Cam5AxisTest, InverseKinematics_AC_TableTable) { + // 创建一个简单的 5 轴刀路 + AxisToolpath5 tp; + tp.name = "Test"; + tp.points.emplace_back(Point3D(10, 0, 0), Vector3D::UnitZ(), 500); + tp.points.emplace_back(Point3D(0, 10, 0), + Vector3D(0, -std::sin(M_PI/4), std::cos(M_PI/4)), 500); + + MachineConfig cfg; + cfg.machine_type = MachineType::TableTable; + cfg.axis_config = MachineAxisConfig::AC; + + auto machine_tp = inverse_kinematics(tp, cfg); + EXPECT_FALSE(machine_tp.empty()); + // 结果应包含机床坐标 + for (const auto& cl : machine_tp.points) { + // tool_axis 存储 A/B/C 度数 + EXPECT_LE(std::abs(cl.tool_axis.x()), 180.0); + EXPECT_LE(std::abs(cl.tool_axis.y()), 180.0); + } +} + +TEST(Cam5AxisTest, InverseKinematics_EmptyInput) { + AxisToolpath5 tp; + tp.name = "Empty"; + MachineConfig cfg; + + auto result = inverse_kinematics(tp, cfg); + EXPECT_TRUE(result.empty()); +} + +// =========================================================================== +// 碰撞检测 check_collision 测试 +// =========================================================================== + +TEST(Cam5AxisTest, CheckCollision_ClearPath) { + auto box = make_test_box2(); + Tool tool; + tool.diameter = 6.0; + tool.length = 20.0; + tool.overall_length = 60.0; + + // 刀尖在模型上方远处，不应碰撞 + Point3D pos(0, 0, 80.0); + Vector3D axis = Vector3D::UnitZ(); + + bool collides = check_collision(pos, axis, box, tool); + EXPECT_FALSE(collides); +} + +TEST(Cam5AxisTest, CheckCollision_InsideModel) { + auto box = make_test_box2(); + Tool tool; + tool.diameter = 6.0; + tool.length = 20.0; + tool.overall_length = 60.0; + + // 刀尖在模型内部 + Point3D pos(0, 0, 15.0); + Vector3D axis = Vector3D::UnitZ(); + + bool collides = check_collision(pos, axis, box, tool); + // 可能检测到碰撞（取决于 mesh tessellation 精度） + // 不做严格断言，只验证函数可调用不崩溃 + SUCCEED(); +} + +// =========================================================================== +// optimize_tool_axis 刀轴优化测试 +// =========================================================================== + +TEST(Cam5AxisTest, OptimizeToolAxis_EmptyCandidates) { + std::vector candidates; + auto result = optimize_tool_axis(candidates, MachineConfig{}); + // 应返回 Z 轴 + EXPECT_NEAR(result.dot(Vector3D::UnitZ()), 1.0, 1e-6); +} + +TEST(Cam5AxisTest, OptimizeToolAxis_PrefersZAxis) { + std::vector candidates = { + Vector3D(0, 0, 1), // Z: 得分最高 + Vector3D(0.5, 0, 0.866), // 30° off + Vector3D(0, -1, 0), // 90° off (poor) + }; + auto result = optimize_tool_axis(candidates, MachineConfig{}); + EXPECT_NEAR(result.dot(Vector3D::UnitZ()), 1.0, 1e-6); +} + +TEST(Cam5AxisTest, OptimizeToolAxis_FiltersOutOfRange) { + MachineConfig cfg; + cfg.a_min = -10.0; + cfg.a_max = 10.0; + cfg.axis_config = MachineAxisConfig::AC; + + // 一个超出范围的方向和一个可行方向 + std::vector candidates = { + Vector3D(0, -1, 0), // 90° tilt → 超出范围 + Vector3D(0, 0, 1), // Z 轴 → 可行 + }; + auto result = optimize_tool_axis(candidates, cfg); + EXPECT_NEAR(result.dot(Vector3D::UnitZ()), 1.0, 1e-6); +} + +// =========================================================================== +// ToolAxisStrategy 枚举完整性测试 +// =========================================================================== + +TEST(Cam5AxisTest, ToolAxisStrategy_AllValues) { + // 确保所有枚举值可被遍历 + EXPECT_EQ(static_cast(ToolAxisStrategy::FixedAngle), 0); + EXPECT_EQ(static_cast(ToolAxisStrategy::LeadLag), 1); + EXPECT_EQ(static_cast(ToolAxisStrategy::Tilted), 2); + EXPECT_EQ(static_cast(ToolAxisStrategy::NormalToSurface), 3); +} + +// =========================================================================== +// MachineConfig 配置完整性测试 +// =========================================================================== + +TEST(Cam5AxisTest, MachineConfig_AllAxisConfigs) { + for (auto ac : {MachineAxisConfig::AC, MachineAxisConfig::BC, MachineAxisConfig::AB}) { + MachineConfig cfg; + cfg.axis_config = ac; + // Z 轴方向应为所有配置都可达 + auto result = cfg.axis_to_angles(Vector3D::UnitZ()); + EXPECT_TRUE(result.has_value()) + << "Axis config failed for Z-up direction"; + } +} + +TEST(Cam5AxisTest, MachineConfig_UnknownConfig) { + // AB 配置：X 轴方向 + MachineConfig cfg; + cfg.axis_config = MachineAxisConfig::AB; + auto result = cfg.axis_to_angles(Vector3D::UnitX()); + ASSERT_TRUE(result.has_value()); + auto [a, b, c] = result.value(); + // AB 配置没有 C 轴 + EXPECT_DOUBLE_EQ(c, 0.0); +} + +// =========================================================================== +// 参数默认值测试 +// =========================================================================== + +TEST(Cam5AxisTest, SwarfParams_Defaults) { + SwarfParams p; + EXPECT_DOUBLE_EQ(p.step_down, 1.0); + EXPECT_DOUBLE_EQ(p.step_over, 0.5); + EXPECT_DOUBLE_EQ(p.feed_rate, 600.0); + EXPECT_DOUBLE_EQ(p.tilt_angle, 0.0); +} + +TEST(Cam5AxisTest, MultiAxisRoughingParams_Defaults) { + MultiAxisRoughingParams p; + EXPECT_DOUBLE_EQ(p.step_down, 2.0); + EXPECT_DOUBLE_EQ(p.stock_to_leave, 1.0); + EXPECT_DOUBLE_EQ(p.tilt_angle, 15.0); +} + +TEST(Cam5AxisTest, MultiAxisFinishingParams_Defaults) { + MultiAxisFinishingParams p; + EXPECT_DOUBLE_EQ(p.step_over, 0.3); + EXPECT_EQ(p.axis_strategy, ToolAxisStrategy::NormalToSurface); + EXPECT_TRUE(p.collision_check); + EXPECT_DOUBLE_EQ(p.shank_clearance, 2.0); + EXPECT_DOUBLE_EQ(p.holder_clearance, 5.0); +} diff --git a/tests/mesh/CMakeLists.txt b/tests/mesh/CMakeLists.txt index 1dfe472..3bc1266 100644 --- a/tests/mesh/CMakeLists.txt +++ b/tests/mesh/CMakeLists.txt @@ -5,3 +5,5 @@ add_vde_test(test_smooth) add_vde_test(test_delaunay_3d) add_vde_test(test_mesh_lod) add_vde_test(test_parallel_mc) +add_vde_test(test_reverse_engineering) +add_vde_test(test_fea_mesh) diff --git a/tests/mesh/test_fea_mesh.cpp b/tests/mesh/test_fea_mesh.cpp new file mode 100644 index 0000000..339c441 --- /dev/null +++ b/tests/mesh/test_fea_mesh.cpp @@ -0,0 +1,356 @@ +#include +#include "vde/mesh/fea_mesh.h" +#include "vde/brep/modeling.h" +#include "vde/mesh/mesh_quality.h" +#include +#include +#include + +using namespace vde::mesh; +using namespace vde::brep; +using namespace vde::core; + +// ═══════════════════════════════════════════════════════════ +// 1. FEAMesh 数据结构测试 +// ═══════════════════════════════════════════════════════════ + +TEST(FEAMeshTest, DefaultConstruction) { + FEAMesh mesh; + EXPECT_EQ(mesh.num_vertices(), 0u); + EXPECT_EQ(mesh.num_elements(), 0u); + EXPECT_EQ(mesh.num_boundary_faces(), 0u); + EXPECT_EQ(mesh.element_type, FEAElementType::Tet4); + EXPECT_EQ(mesh.npe(), 4); +} + +TEST(FEAMeshTest, ElementTypeNPE) { + EXPECT_EQ(nodes_per_element(FEAElementType::Tet4), 4); + EXPECT_EQ(nodes_per_element(FEAElementType::Tet10), 10); + EXPECT_EQ(nodes_per_element(FEAElementType::Hex8), 8); + EXPECT_EQ(nodes_per_element(FEAElementType::Hex20), 20); + EXPECT_EQ(nodes_per_element(FEAElementType::Wedge6), 6); + EXPECT_EQ(nodes_per_element(FEAElementType::Wedge15), 15); +} + +TEST(FEAMeshTest, ElementTypeName) { + EXPECT_STREQ(element_type_name(FEAElementType::Tet4), "Tet4"); + EXPECT_STREQ(element_type_name(FEAElementType::Hex8), "Hex8"); + EXPECT_STREQ(element_type_name(FEAElementType::Wedge6), "Wedge6"); +} + +// ═══════════════════════════════════════════════════════════ +// 2. tetrahedral_mesh 测试 +// ═══════════════════════════════════════════════════════════ + +TEST(TetrahedralMeshTest, BoxMesh_NonEmpty) { + auto box = make_box(1.0, 1.0, 1.0); + TetMeshParams params; + params.max_size = 0.3; + params.quality_iterations = 1; + + auto mesh = tetrahedral_mesh(box, params); + EXPECT_GT(mesh.num_vertices(), 0u); + EXPECT_GT(mesh.num_elements(), 0u); + EXPECT_EQ(mesh.element_type, FEAElementType::Tet4); + // 四面体网格应该有边界面 + EXPECT_GT(mesh.num_boundary_faces(), 0u); +} + +TEST(TetrahedralMeshTest, BoxMesh_ValidConnectivity) { + auto box = make_box(1.0, 1.0, 1.0); + TetMeshParams params; + params.max_size = 0.5; + + auto mesh = tetrahedral_mesh(box, params); + // 所有单元的顶点索引应在合法范围内 + for (size_t ei = 0; ei < mesh.num_elements(); ++ei) { + auto& e = mesh.elements[ei]; + EXPECT_EQ(e.size(), 4u); + for (int vi : e) { + EXPECT_GE(vi, 0); + EXPECT_LT(static_cast(vi), mesh.num_vertices()); + } + } +} + +TEST(TetrahedralMeshTest, SphereMesh_NonTrivial) { + auto sphere = make_sphere(1.0); + TetMeshParams params; + params.max_size = 0.5; + params.quality_iterations = 1; + + auto mesh = tetrahedral_mesh(sphere, params); + EXPECT_GT(mesh.num_vertices(), 0u); + EXPECT_GT(mesh.num_elements(), 0u); +} + +// ═══════════════════════════════════════════════════════════ +// 3. boundary_layer_mesh 测试 +// ═══════════════════════════════════════════════════════════ + +TEST(BoundaryLayerTest, BoxBoundaryLayer_Prisms) { + auto box = make_box(1.0, 1.0, 1.0); + BLPParams params; + params.first_cell_height = 0.01; + params.growth_rate = 1.2; + params.num_layers = 3; + + auto mesh = boundary_layer_mesh(box, params); + EXPECT_EQ(mesh.element_type, FEAElementType::Wedge6); + EXPECT_GT(mesh.num_elements(), 0u); + EXPECT_GT(mesh.num_vertices(), 0u); + // 顶点数 = 表面顶点 × (layers+1) + EXPECT_GT(mesh.num_vertices(), 0u); +} + +TEST(BoundaryLayerTest, BoxBoundaryLayer_ValidConnectivity) { + auto box = make_box(1.0, 1.0, 1.0); + BLPParams params; + params.num_layers = 2; + + auto mesh = boundary_layer_mesh(box, params); + for (size_t ei = 0; ei < mesh.num_elements(); ++ei) { + auto& e = mesh.elements[ei]; + EXPECT_EQ(e.size(), 6u); // Wedge6 + for (int vi : e) { + EXPECT_GE(vi, 0); + EXPECT_LT(static_cast(vi), mesh.num_vertices()); + } + } +} + +// ═══════════════════════════════════════════════════════════ +// 4. hexahedral_mesh 测试 +// ═══════════════════════════════════════════════════════════ + +TEST(HexahedralMeshTest, BoxSweep_HexMesh) { + auto box = make_box(1.0, 1.0, 1.0); + HexMeshParams params; + params.sweep_layers = 2; + + auto mesh = hexahedral_mesh(box, params); + EXPECT_EQ(mesh.element_type, FEAElementType::Hex8); + EXPECT_GE(mesh.num_elements(), 0u); + EXPECT_GT(mesh.num_vertices(), 0u); +} + +TEST(HexahedralMeshTest, BoxSweep_ValidConnectivity) { + auto box = make_box(1.0, 1.0, 1.0); + HexMeshParams params; + params.sweep_layers = 2; + + auto mesh = hexahedral_mesh(box, params); + for (size_t ei = 0; ei < mesh.num_elements(); ++ei) { + auto& e = mesh.elements[ei]; + EXPECT_EQ(e.size(), 8u); // Hex8 + for (int vi : e) { + EXPECT_GE(vi, 0); + EXPECT_LT(static_cast(vi), mesh.num_vertices()); + } + } +} + +// ═══════════════════════════════════════════════════════════ +// 5. 单元素质量指标测试 +// ═══════════════════════════════════════════════════════════ + +TEST(ElementQualityTest, Tet4_Regular_GoodQuality) { + // 正四面体 (边长为 sqrt(2) 的四个点) + std::vector verts = { + {1, 1, 1}, + {1, -1, -1}, + {-1, 1, -1}, + {-1, -1, 1} + }; + double sj = element_scaled_jacobian(verts, FEAElementType::Tet4); + EXPECT_GT(sj, 0.5); + + double skew = element_skewness(verts, FEAElementType::Tet4); + EXPECT_LT(skew, 0.5); + + double ortho = element_orthogonality(verts, FEAElementType::Tet4); + EXPECT_GT(ortho, 0.0); + EXPECT_LE(ortho, 1.0); +} + +TEST(ElementQualityTest, Tet4_Degenerate_ZeroJacobian) { + // 退化四面体（四点共面） + std::vector verts = { + {0, 0, 0}, + {1, 0, 0}, + {0, 1, 0}, + {0.5, 0.5, 0} // 在同一平面上 + }; + double sj = element_scaled_jacobian(verts, FEAElementType::Tet4); + EXPECT_NEAR(sj, 0.0, 1e-6); +} + +TEST(ElementQualityTest, AspectRatio_IsotropicElement) { + // 边长为 1 的四面体 + std::vector verts = { + {0, 0, 0}, + {1, 0, 0}, + {0, 1, 0}, + {0, 0, 1} + }; + double ar = element_aspect_ratio(verts, FEAElementType::Tet4); + EXPECT_GE(ar, 1.0); +} + +TEST(ElementQualityTest, Hex8_QualityFinite) { + std::vector verts = { + {0,0,0},{1,0,0},{1,1,0},{0,1,0}, + {0,0,1},{1,0,1},{1,1,1},{0,1,1} + }; + double sj = element_scaled_jacobian(verts, FEAElementType::Hex8); + EXPECT_GT(sj, 0.0); + EXPECT_LE(sj, 1.0); + + double skew = element_skewness(verts, FEAElementType::Hex8); + EXPECT_LE(skew, 1.0); +} + +// ═══════════════════════════════════════════════════════════ +// 6. fea_quality_report 测试 +// ═══════════════════════════════════════════════════════════ + +TEST(FEAQualityReportTest, EmptyMesh_AllZero) { + FEAMesh mesh; + auto report = fea_quality_report(mesh); + EXPECT_EQ(report.total_elements, 0u); + EXPECT_EQ(report.degenerate_elements, 0u); +} + +TEST(FEAQualityReportTest, TetrahedralMesh_ReportValid) { + auto box = make_box(1.0, 1.0, 1.0); + TetMeshParams params; + params.max_size = 0.4; + auto mesh = tetrahedral_mesh(box, params); + + auto report = fea_quality_report(mesh); + EXPECT_EQ(report.total_elements, mesh.num_elements()); + EXPECT_GT(report.total_elements, 0u); + // 质量报告应有有效值 + EXPECT_GE(report.avg_jacobian, 0.0); + EXPECT_LE(report.avg_jacobian, 1.0); + EXPECT_GE(report.avg_skewness, 0.0); + EXPECT_LE(report.avg_skewness, 1.0); +} + +// ═══════════════════════════════════════════════════════════ +// 7. adaptive_refinement 测试 +// ═══════════════════════════════════════════════════════════ + +TEST(AdaptiveRefinementTest, NoRefinement_ReturnsSame) { + auto box = make_box(1.0, 1.0, 1.0); + auto mesh = tetrahedral_mesh(box, TetMeshParams{}); + size_t original_elems = mesh.num_elements(); + + // 误差为 0 → 不细化 + auto zero_estimator = [](int, const FEAMesh&) -> double { return 0.0; }; + auto refined = adaptive_refinement(mesh, zero_estimator); + + // 不细化时单元数应不变（但边中点缓存可能导致微小差异） + // 只验证不崩溃且仍有单元 + EXPECT_GT(refined.num_elements(), 0u); +} + +TEST(AdaptiveRefinementTest, HighError_Refines) { + auto box = make_box(1.0, 1.0, 1.0); + TetMeshParams params; + params.max_size = 0.5; + auto mesh = tetrahedral_mesh(box, params); + size_t original_elems = mesh.num_elements(); + + // 所有单元高误差 → 细化 + auto high_estimator = [](int, const FEAMesh&) -> double { return 1.0; }; + RefinementParams rp; + rp.error_threshold = 0.1; + auto refined = adaptive_refinement(mesh, high_estimator, rp); + + // 细化后单元数应增加 + EXPECT_GT(refined.num_elements(), original_elems); +} + +// ═══════════════════════════════════════════════════════════ +// 8. CAE 导出测试 +// ═══════════════════════════════════════════════════════════ + +TEST(CAEExportTest, AbaqusExport_FileCreated) { + auto box = make_box(1.0, 1.0, 1.0); + auto mesh = tetrahedral_mesh(box, TetMeshParams{}); + + std::string path = "/tmp/test_abaqus.inp"; + bool ok = export_abaqus(mesh, path); + EXPECT_TRUE(ok); + + // 检查文件存在且非空 + std::ifstream f(path); + EXPECT_TRUE(f.good()); + std::string content((std::istreambuf_iterator(f)), + std::istreambuf_iterator()); + EXPECT_GT(content.size(), 0u); + // 应包含关键关键字 + EXPECT_NE(content.find("*NODE"), std::string::npos); + EXPECT_NE(content.find("*ELEMENT"), std::string::npos); + std::remove(path.c_str()); +} + +TEST(CAEExportTest, AnsysExport_FileCreated) { + auto box = make_box(1.0, 1.0, 1.0); + auto mesh = tetrahedral_mesh(box, TetMeshParams{}); + + std::string path = "/tmp/test_ansys.cdb"; + bool ok = export_ansys(mesh, path); + EXPECT_TRUE(ok); + + std::ifstream f(path); + EXPECT_TRUE(f.good()); + std::string content((std::istreambuf_iterator(f)), + std::istreambuf_iterator()); + EXPECT_GT(content.size(), 0u); + EXPECT_NE(content.find("NBLOCK"), std::string::npos); + EXPECT_NE(content.find("EBLOCK"), std::string::npos); + std::remove(path.c_str()); +} + +TEST(CAEExportTest, NastranExport_FileCreated) { + auto box = make_box(1.0, 1.0, 1.0); + auto mesh = tetrahedral_mesh(box, TetMeshParams{}); + + std::string path = "/tmp/test_nastran.bdf"; + bool ok = export_nastran(mesh, path); + EXPECT_TRUE(ok); + + std::ifstream f(path); + EXPECT_TRUE(f.good()); + std::string content((std::istreambuf_iterator(f)), + std::istreambuf_iterator()); + EXPECT_GT(content.size(), 0u); + EXPECT_NE(content.find("GRID"), std::string::npos); + EXPECT_NE(content.find("CTETRA"), std::string::npos); + std::remove(path.c_str()); +} + +TEST(CAEExportTest, Export_InvalidPath) { + FEAMesh mesh; + bool ok = export_abaqus(mesh, "/nonexistent_dir/should_fail.inp"); + EXPECT_FALSE(ok); +} + +TEST(CAEExportTest, HexMeshNastranExport) { + auto box = make_box(2.0, 1.0, 1.0); + HexMeshParams params; + params.sweep_layers = 2; + auto mesh = hexahedral_mesh(box, params); + + std::string path = "/tmp/test_hex_nastran.bdf"; + bool ok = export_nastran(mesh, path); + EXPECT_TRUE(ok); + + std::ifstream f(path); + std::string content((std::istreambuf_iterator(f)), + std::istreambuf_iterator()); + EXPECT_NE(content.find("CHEXA"), std::string::npos); + std::remove(path.c_str()); +} diff --git a/tests/mesh/test_reverse_engineering.cpp b/tests/mesh/test_reverse_engineering.cpp new file mode 100644 index 0000000..08c34fa --- /dev/null +++ b/tests/mesh/test_reverse_engineering.cpp @@ -0,0 +1,268 @@ +#include +#include "vde/mesh/reverse_engineering.h" +#include "vde/mesh/point_cloud.h" +#include "vde/mesh/halfedge_mesh.h" +#include "vde/curves/nurbs_surface.h" + +using namespace vde::mesh; +using namespace vde::curves; + +// ═══════════════════════════════════════════════════════════ +// Helper: generate a simple point cloud (e.g., a noisy plane) +// ═══════════════════════════════════════════════════════════ + +static PointCloud make_plane_cloud(int n = 100) { + std::vector pts; + for (int i = 0; i < n; ++i) { + double u = static_cast(i % 10) / 10.0; + double v = static_cast(i / 10) / 10.0; + pts.emplace_back(u, v, 0.01 * std::sin(u * 10) * std::cos(v * 10)); + } + return PointCloud(pts); +} + +static PointCloud make_sphere_cloud(int n = 200) { + std::vector pts; + for (int i = 0; i < n; ++i) { + double theta = 2.0 * M_PI * static_cast(i) / n; + double phi = M_PI * static_cast((i * 7) % n) / n; + double x = std::sin(phi) * std::cos(theta); + double y = std::sin(phi) * std::sin(theta); + double z = std::cos(phi); + pts.emplace_back(x, y, z); + } + return PointCloud(pts); +} + +// ═══════════════════════════════════════════════════════════ +// PointCloud 基础测试 +// ═══════════════════════════════════════════════════════════ + +TEST(PointCloudTest, Construction) { + PointCloud empty; + EXPECT_EQ(empty.size(), 0u); + EXPECT_FALSE(empty.has_normals()); + EXPECT_FALSE(empty.has_colors()); + + PointCloud cloud(std::vector{{0,0,0}, {1,0,0}, {0,1,0}}); + EXPECT_EQ(cloud.size(), 3u); +} + +TEST(PointCloudTest, Bounds) { + PointCloud cloud(std::vector{{0,0,0}, {2,0,0}, {0,3,0}, {2,3,0}}); + auto bb = cloud.bounds(); + EXPECT_NEAR(bb.min().x(), 0.0, 1e-9); + EXPECT_NEAR(bb.min().y(), 0.0, 1e-9); + EXPECT_NEAR(bb.max().x(), 2.0, 1e-9); + EXPECT_NEAR(bb.max().y(), 3.0, 1e-9); +} + +TEST(PointCloudTest, VoxelDownsample) { + std::vector pts; + for (int i = 0; i < 50; ++i) pts.emplace_back(i * 0.01, 0, 0); + PointCloud cloud(pts); + auto down = cloud.voxel_downsample(0.1); + EXPECT_GT(down.size(), 0u); + EXPECT_LT(down.size(), pts.size()); +} + +TEST(PointCloudTest, StatisticalOutlierRemoval) { + std::vector pts; + // Dense cluster + for (int i = 0; i < 100; ++i) + pts.emplace_back(0.01 * i, 0, 0); + // Outlier far away + pts.emplace_back(100, 100, 100); + PointCloud cloud(pts); + auto clean = cloud.statistical_outlier_removal(6, 1.0); + EXPECT_LT(clean.size(), pts.size()); + // The outlier should be removed + for (const auto& p : clean.vertices) { + EXPECT_LT((p - Point3D(100, 100, 100)).norm(), 50.0); + } +} + +TEST(PointCloudTest, NormalEstimation) { + auto cloud = make_plane_cloud(100); + auto with_normals = cloud.normal_estimation(10); + EXPECT_EQ(with_normals.size(), cloud.size()); + EXPECT_TRUE(with_normals.has_normals()); + // Plane normals should be mostly Z-aligned + for (size_t i = 0; i < with_normals.size(); ++i) { + double dot_z = std::abs(with_normals.normals[i].z()); + EXPECT_GT(dot_z, 0.7); // Should be mostly pointing in Z + } +} + +TEST(PointCloudTest, SmoothingFilter) { + std::vector pts; + for (int i = 0; i < 50; ++i) + pts.emplace_back(i * 0.02, 0, 0.1 * std::sin(i * 0.5)); + PointCloud cloud(pts); + auto smoothed = cloud.smoothing_filter(0.5, 0.3, 3); + EXPECT_EQ(smoothed.size(), cloud.size()); +} + +TEST(PointCloudTest, CurvatureAwareSampling) { + auto cloud = make_sphere_cloud(200); + auto indices = cloud.curvature_aware_sampling(50); + EXPECT_LE(indices.size(), 50u); + EXPECT_GT(indices.size(), 0u); +} + +// ═══════════════════════════════════════════════════════════ +// Plane fitting +// ═══════════════════════════════════════════════════════════ + +TEST(ReverseEngineeringTest, FitPlane_XY_Plane) { + std::vector pts; + for (int i = 0; i < 100; ++i) + pts.emplace_back(static_cast(i % 10), static_cast(i / 10), 0.0); + + auto result = fit_plane(pts); + EXPECT_TRUE(result.valid); + EXPECT_NEAR(result.rms_error, 0.0, 1e-6); + EXPECT_NEAR(std::abs(result.normal.z()), 1.0, 1e-6); +} + +TEST(ReverseEngineeringTest, FitPlane_Slanted) { + std::vector pts; + // z = x + y + for (int i = 0; i < 100; ++i) { + double x = static_cast(i % 10); + double y = static_cast(i / 10); + pts.emplace_back(x, y, x + y); + } + + auto result = fit_plane(pts); + EXPECT_TRUE(result.valid); + EXPECT_NEAR(result.rms_error, 0.0, 1e-6); + // Normal should be perpendicular to (1,1,-1) → normalize + double dot = result.normal.x() + result.normal.y() - result.normal.z(); + EXPECT_NEAR(std::abs(dot), std::sqrt(3.0), 1e-6); +} + +TEST(ReverseEngineeringTest, FitPlane_InsufficientPoints) { + std::vector pts = {{0,0,0}, {1,0,0}}; + auto result = fit_plane(pts); + EXPECT_FALSE(result.valid); +} + +// ═══════════════════════════════════════════════════════════ +// Patch / NURBS fitting +// ═══════════════════════════════════════════════════════════ + +TEST(ReverseEngineeringTest, FitPatch_PlaneCloud) { + auto cloud = make_plane_cloud(100); + auto result = fit_patch(cloud.vertices, 2, 2, 5, 5); + EXPECT_TRUE(result.valid); + EXPECT_GT(result.rms_error, 0.0); + // Surface should be evaluable + auto pt = result.surface.evaluate(0.5, 0.5); + EXPECT_TRUE(std::isfinite(pt.x())); +} + +TEST(ReverseEngineeringTest, PointCloudToMesh) { + auto cloud = make_plane_cloud(200); + PointCloudToMeshParams p; + p.knn = 8; + p.poisson.max_depth = 5; + auto mesh = point_cloud_to_mesh(cloud, p); + EXPECT_GT(mesh.num_vertices(), 0u); + EXPECT_GT(mesh.num_faces(), 0u); +} + +TEST(ReverseEngineeringTest, MeshToNurbsSurface) { + // Build simple mesh from plane cloud + auto cloud = make_plane_cloud(100); + PointCloudToMeshParams pm; + pm.knn = 8; + pm.poisson.max_depth = 5; + auto mesh = point_cloud_to_mesh(cloud, pm); + + ASSERT_GT(mesh.num_vertices(), 0u); + + MeshToNURBSParams mp; + mp.nurbs.num_cp_u = 5; + mp.nurbs.num_cp_v = 5; + mp.nurbs.degree_u = 2; + mp.nurbs.degree_v = 2; + mp.nurbs.max_iterations = 10; + auto surf = mesh_to_nurbs_surface(mesh, mp); + // Surface should be valid + auto pt = surf.evaluate(0.5, 0.5); + EXPECT_TRUE(std::isfinite(pt.x())); +} + +// ═══════════════════════════════════════════════════════════ +// Curvature-aware sampling (free function) +// ═══════════════════════════════════════════════════════════ + +TEST(ReverseEngineeringTest, CurvatureAwareSampling_Basic) { + auto cloud = make_sphere_cloud(200); + auto indices = curvature_aware_sampling(cloud.vertices, 50); + EXPECT_LE(indices.size(), 50u); + EXPECT_GT(indices.size(), 0u); + // Indices should be within range + for (auto idx : indices) EXPECT_LT(idx, cloud.vertices.size()); +} + +TEST(ReverseEngineeringTest, CurvatureAwareSampling_AllPoints) { + auto cloud = make_plane_cloud(30); + auto indices = curvature_aware_sampling(cloud.vertices, 100); + EXPECT_EQ(indices.size(), cloud.vertices.size()); +} + +// ═══════════════════════════════════════════════════════════ +// Hole filling +// ═══════════════════════════════════════════════════════════ + +TEST(ReverseEngineeringTest, HoleFilling_NoHoles) { + // Create a closed mesh (tetrahedron-like) + HalfedgeMesh mesh; + mesh.build_from_triangles( + {{0,0,0}, {1,0,0}, {0,1,0}, {0,0,1}}, + {{0,1,2}, {0,1,3}, {1,2,3}, {0,2,3}} + ); + auto filled = hole_filling(mesh); + // Should be same size or larger + EXPECT_GE(filled.num_faces(), mesh.num_faces()); +} + +// ═══════════════════════════════════════════════════════════ +// Smoothing filter (free function) +// ═══════════════════════════════════════════════════════════ + +TEST(ReverseEngineeringTest, SmoothingFilter_Basic) { + auto cloud = make_plane_cloud(100); + auto smoothed = smoothing_filter(cloud.vertices, 0.5, 3); + EXPECT_EQ(smoothed.size(), cloud.vertices.size()); +} + +TEST(ReverseEngineeringTest, SmoothingFilter_EmptyPoints) { + std::vector empty; + auto smoothed = smoothing_filter(empty, 0.5, 3); + EXPECT_TRUE(smoothed.empty()); +} + +// ═══════════════════════════════════════════════════════════ +// Estimate normals (free function) +// ═══════════════════════════════════════════════════════════ + +TEST(ReverseEngineeringTest, EstimateNormals_Plane) { + std::vector pts; + for (int i = 0; i < 50; ++i) + pts.emplace_back(i * 0.02, 0.0, 0.0); + pts.emplace_back(0.0, 0.02, 0.0); + for (int i = 0; i < 50; ++i) + pts.emplace_back(0.0, i * 0.02, 0.0); + + auto normals = estimate_normals(pts, 10); + EXPECT_EQ(normals.size(), pts.size()); + // Most normals should be Z-aligned + int z_count = 0; + for (const auto& n : normals) { + if (std::abs(n.z()) > 0.7) ++z_count; + } + EXPECT_GT(z_count, static_cast(pts.size() * 0.8)); +}