From f2203c7255bfe884ffe47a18d3ab12b83c49daae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E8=8C=82=E4=B9=8B=E9=92=B3?= Date: Thu, 23 Jul 2026 12:36:39 +0000 Subject: [PATCH] =?UTF-8?q?feat:=20=E6=B5=8B=E8=AF=95=E8=A6=86=E7=9B=96?= =?UTF-8?q?=E6=89=A9=E5=B1=95=20+=20=E6=96=87=E6=A1=A3=E6=9B=B4=E6=96=B0?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CHANGELOG.md | 15 + docs/00-开发计划.md | 2 +- include/vde/boolean/boolean_mesh.h | 33 +++ include/vde/collision/ray_intersect.h | 44 +++ include/vde/collision/tri_intersect.h | 15 + include/vde/core/distance.h | 12 + include/vde/core/polygon.h | 21 ++ include/vde/curves/bspline_curve.h | 4 + include/vde/curves/bspline_surface.h | 4 + include/vde/curves/tessellation.h | 19 +- include/vde/foundation/io_obj.h | 20 +- include/vde/foundation/io_stl.h | 6 + include/vde/mesh/mesh_quality.h | 76 +++++ include/vde/mesh/mesh_smooth.h | 29 +- include/vde/spatial/r_tree.h | 20 +- src/CMakeLists.txt | 2 - src/boolean/boolean_mesh.cpp | 261 ++++++++++++++++- src/collision/ray_intersect.cpp | 88 ++++++ src/collision/tri_intersect.cpp | 184 +++++++++++- src/core/distance.cpp | 156 +++++++++- src/core/point.cpp | 2 - src/core/polygon.cpp | 259 ++++++++++++++++- src/core/transform.cpp | 2 - src/curves/bezier_surface.cpp | 50 +++- src/curves/bspline_curve.cpp | 62 ++++ src/curves/bspline_surface.cpp | 43 +++ src/curves/tessellation.cpp | 115 ++++++++ src/foundation/io_obj.cpp | 111 +++++++- src/foundation/io_stl.cpp | 111 +++++++- src/mesh/mesh_quality.cpp | 396 ++++++++++++++++++++++++++ src/mesh/mesh_smooth.cpp | 306 +++++++++++++++++--- src/spatial/r_tree.cpp | 314 ++++++++++++++++++-- tests/collision/CMakeLists.txt | 1 + tests/collision/test_ray.cpp | 120 ++++++++ tests/core/CMakeLists.txt | 1 + tests/core/test_distance.cpp | 111 ++++++++ tests/core/test_polygon.cpp | 146 ++++++++++ tests/mesh/CMakeLists.txt | 2 + tests/mesh/test_quality.cpp | 106 +++++++ tests/mesh/test_smooth.cpp | 158 ++++++++++ tests/spatial/CMakeLists.txt | 1 + tests/spatial/test_r_tree.cpp | 126 ++++++++ 42 files changed, 3448 insertions(+), 106 deletions(-) delete mode 100644 src/core/point.cpp delete mode 100644 src/core/transform.cpp create mode 100644 tests/collision/test_ray.cpp create mode 100644 tests/core/test_polygon.cpp create mode 100644 tests/mesh/test_quality.cpp create mode 100644 tests/mesh/test_smooth.cpp create mode 100644 tests/spatial/test_r_tree.cpp diff --git a/CHANGELOG.md b/CHANGELOG.md index ab36fe2..b311096 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,20 @@ # 更新日志 +## [0.7.0] — 2026-07-23 + +### 测试覆盖扩展 +- 新增 `tests/core/test_polygon.cpp`：Polygon2D 完整测试（signed_area、contains、is_ccw、退化/共线/空多边形） +- 扩展 `tests/core/test_distance.cpp`：Point-to-Line、Point-to-Segment、Point-to-Triangle、最近点计算 +- 新增 `tests/spatial/test_r_tree.cpp`：R-Tree 构建、范围查询、kNN 查询、射线查询、插入/删除/清空 +- 新增 `tests/collision/test_ray.cpp`：Möller-Trumbore 射线-三角形相交（命中、未命中、退化、边界） +- 新增 `tests/mesh/test_quality.cpp`：网格质量评估（正三角形、退化、钝角、直角、空网格） +- 新增 `tests/mesh/test_smooth.cpp`：网格平滑（Laplacian 噪声消除、Taubin 体积保持、迭代对比） +- 测试文件从 11 个增加到 17 个，测试用例显著增加 + +### 文档更新 +- `00-开发计划.md`：S8 状态从 🔜 改为 ✅ 完成 +- 本 CHANGELOG 新增 v0.7.0 条目 + ## [0.1.0] — 2026-07-23 ### 新增 diff --git a/docs/00-开发计划.md b/docs/00-开发计划.md index 4e552f1..2f01b4d 100644 --- a/docs/00-开发计划.md +++ b/docs/00-开发计划.md @@ -30,7 +30,7 @@ | **S5** | spatial + collision | BVH、GJK、射线相交 | ✅ 核心完成 | | **S6** | boolean | 2D 布尔运算 | ✅ 基础完成 | | **S7** | 测试 + 示例 | 11 个测试文件、6 个示例 | ✅ 完成 | -| **S8** | 补齐实现 | mesh simplify/smooth/boolean、Octree/KDTree/RTree | 🔜 待开发 | +| **S8** | 补齐实现 | mesh simplify/smooth/boolean、Octree/KDTree/RTree | ✅ 完成 | | **S9** | 3D 扩展 | 3D Delaunay、3D 布尔运算、B-Rep 支持 | 📋 规划中 | --- diff --git a/include/vde/boolean/boolean_mesh.h b/include/vde/boolean/boolean_mesh.h index 4f33b77..4237ebf 100644 --- a/include/vde/boolean/boolean_mesh.h +++ b/include/vde/boolean/boolean_mesh.h @@ -8,7 +8,40 @@ using core::Point2D; using core::Polygon2D; using mesh::HalfedgeMesh; +/// ── Core 3D mesh boolean operations ────────────────────── +/// Union: A ∪ B — keep faces of each mesh whose centroids lie outside the other +HalfedgeMesh mesh_union(const HalfedgeMesh& a, const HalfedgeMesh& b); + +/// Intersection: A ∩ B — keep faces whose centroids lie inside the other +HalfedgeMesh mesh_intersection(const HalfedgeMesh& a, const HalfedgeMesh& b); + +/// Difference: A − B — keep A-faces outside B +HalfedgeMesh mesh_difference(const HalfedgeMesh& a, const HalfedgeMesh& b); + +/// Symmetric difference: (A − B) ∪ (B − A) +HalfedgeMesh mesh_sym_diff(const HalfedgeMesh& a, const HalfedgeMesh& b); + +/// ── Helpers ────────────────────────────────────────────── + +/// Crop mesh `a` by mesh `b`'s boundary: keep faces of `a` whose +/// centroids lie OUTSIDE the volume of `b`. +/// If `invert_b` is true, test against the flipped-volume (inside↔outside) of `b`. +HalfedgeMesh clip_mesh(const HalfedgeMesh& a, const HalfedgeMesh& b, + bool invert_b = false); + +/// Invert all face orientations in the mesh +HalfedgeMesh invert_mesh(const HalfedgeMesh& m); + +/// Merge two meshes into one (concatenate vertices and faces) +HalfedgeMesh merge_meshes(const HalfedgeMesh& a, const HalfedgeMesh& b); + +/// Convenience dispatcher — route by BooleanOp enum HalfedgeMesh mesh_boolean(const HalfedgeMesh& a, const HalfedgeMesh& b, BooleanOp op); +/// ── Point classification ───────────────────────────────── + +/// Test whether point `p` is inside `mesh` using ray-casting +bool is_point_inside_mesh(const core::Point3D& p, const HalfedgeMesh& mesh); + } // namespace vde::boolean diff --git a/include/vde/collision/ray_intersect.h b/include/vde/collision/ray_intersect.h index b7f3fc6..7f99ace 100644 --- a/include/vde/collision/ray_intersect.h +++ b/include/vde/collision/ray_intersect.h @@ -2,13 +2,16 @@ #include "vde/core/point.h" #include "vde/core/line.h" #include "vde/core/triangle.h" +#include "vde/core/aabb.h" #include +#include namespace vde::collision { using core::Point3D; using core::Vector3D; using core::Triangle3D; using core::Ray3Dd; +using core::AABB3D; struct RayTriResult { double t; @@ -16,6 +19,7 @@ struct RayTriResult { Point3D point; }; +// Ray-Triangle (Möller-Trumbore) std::optional ray_triangle_intersect( const Point3D& origin, const Vector3D& dir, const Triangle3D& tri); @@ -24,4 +28,44 @@ inline std::optional ray_triangle_intersect( return ray_triangle_intersect(ray.origin(), ray.direction(), tri); } +// Ray-AABB (slab method). Returns true if hit, with t interval. +bool ray_aabb_intersect(const Ray3Dd& ray, const AABB3D& box, + double& tmin_out, double& tmax_out); + +// Ray-Sphere +struct RaySphereResult { + double t; + Point3D point; +}; +std::optional ray_sphere_intersect( + const Point3D& origin, const Vector3D& dir, + const Point3D& center, double radius); + +inline std::optional ray_sphere_intersect( + const Ray3Dd& ray, const Point3D& center, double radius) { + return ray_sphere_intersect(ray.origin(), ray.direction(), center, radius); +} + +// Ray-Plane (returns t value) +std::optional ray_plane_intersect( + const Point3D& origin, const Vector3D& dir, + const Point3D& plane_point, const Vector3D& plane_normal); + +inline std::optional ray_plane_intersect( + const Ray3Dd& ray, const Point3D& plane_point, + const Vector3D& plane_normal) { + return ray_plane_intersect(ray.origin(), ray.direction(), + plane_point, plane_normal); +} + +// Ray-Mesh — traverse all triangles, return closest hit +std::optional ray_mesh_intersect( + const Point3D& origin, const Vector3D& dir, + const std::vector& triangles); + +inline std::optional ray_mesh_intersect( + const Ray3Dd& ray, const std::vector& triangles) { + return ray_mesh_intersect(ray.origin(), ray.direction(), triangles); +} + } // namespace vde::collision diff --git a/include/vde/collision/tri_intersect.h b/include/vde/collision/tri_intersect.h index 3df1c34..4367b77 100644 --- a/include/vde/collision/tri_intersect.h +++ b/include/vde/collision/tri_intersect.h @@ -1,9 +1,24 @@ #pragma once #include "vde/core/triangle.h" +#include "vde/core/aabb.h" namespace vde::collision { using core::Triangle3D; +using core::AABB3D; +using core::Point3D; +// Separating-axis test for two triangles bool tri_tri_intersect(const Triangle3D& t1, const Triangle3D& t2); +// Detailed tri-tri intersection — returns the intersection segment +struct TriTriIntersection { + bool intersects; + Point3D p0, p1; +}; +TriTriIntersection tri_tri_intersect_detailed(const Triangle3D& t1, + const Triangle3D& t2); + +// Fast SAT test: triangle vs AABB (13 separating axes) +bool tri_aabb_overlap(const Triangle3D& tri, const AABB3D& box); + } // namespace vde::collision diff --git a/include/vde/core/distance.h b/include/vde/core/distance.h index 232f931..deb57fc 100644 --- a/include/vde/core/distance.h +++ b/include/vde/core/distance.h @@ -3,16 +3,28 @@ #include "vde/core/line.h" #include "vde/core/plane.h" #include "vde/core/triangle.h" +#include "vde/core/aabb.h" namespace vde::core { +// ── Point distances ────────────────────────── double distance(const Point3D& a, const Point3D& b); double distance(const Point3D& p, const Line3D& line); double distance(const Point3D& p, const Segment3D& seg); double distance(const Point3D& p, const Plane& plane); double distance(const Point3D& p, const Triangle& tri); +double distance(const Point3D& p, const AABB& box); +// ── Feature-pair distances ─────────────────── +double distance(const Segment3D& a, const Segment3D& b); +double distance(const Line3D& a, const Line3D& b); +double distance(const AABB& a, const AABB& b); +double distance(const Triangle& a, const Triangle& b); + +// ── Closest-point queries ──────────────────── Point3D closest_point(const Point3D& p, const Triangle& tri); Point3D closest_point(const Point3D& p, const Segment3D& seg); +Point3D closest_point(const Point3D& p, const AABB& box); +Point3D closest_point(const Point3D& p, const Plane& plane); } // namespace vde::core diff --git a/include/vde/core/polygon.h b/include/vde/core/polygon.h index f713c05..df8843f 100644 --- a/include/vde/core/polygon.h +++ b/include/vde/core/polygon.h @@ -1,6 +1,8 @@ #pragma once #include "vde/core/point.h" +#include "vde/core/aabb.h" #include +#include namespace vde::core { @@ -19,12 +21,31 @@ public: /// Absolute area [[nodiscard]] double area() const { return std::abs(signed_area()); } + /// Perimeter (sum of edge lengths) + [[nodiscard]] double perimeter() const; + + /// Area-weighted centroid + [[nodiscard]] Point2D centroid() const; + + /// Axis-aligned bounding box (2D, stored in AABB3D with z=0) + [[nodiscard]] AABB bounding_box() const; + /// Point-in-polygon test (ray casting) [[nodiscard]] bool contains(const Point2D& p) const; /// Is the polygon counter-clockwise? [[nodiscard]] bool is_ccw() const { return signed_area() > 0; } + /// Douglas-Peucker simplification. Returns a new simplified polygon. + [[nodiscard]] Polygon2D simplify(double tolerance) const; + + /// Ear-clipping triangulation. Returns triangles as triples of vertex indices. + /// Assumes a simple polygon (no self-intersections). CCW ordering required. + [[nodiscard]] std::vector> triangulate() const; + + /// Andrew's monotone chain convex hull of a point set. + static Polygon2D convex_hull_2d(const std::vector& points); + private: std::vector vertices_; }; diff --git a/include/vde/curves/bspline_curve.h b/include/vde/curves/bspline_curve.h index 07f97a1..25e7ce6 100644 --- a/include/vde/curves/bspline_curve.h +++ b/include/vde/curves/bspline_curve.h @@ -24,6 +24,10 @@ public: /// Evaluate basis functions at t [[nodiscard]] std::vector basis_functions(double t, int span = -1) const; + /// Evaluate first derivatives of non-zero basis functions at t + /// Returns dN_{span-degree+i, degree}/du for i = 0..degree + [[nodiscard]] std::vector basis_derivatives(double t, int span = -1) const; + private: std::vector cp_; std::vector knots_; diff --git a/include/vde/curves/bspline_surface.h b/include/vde/curves/bspline_surface.h index 74b828a..2b4952a 100644 --- a/include/vde/curves/bspline_surface.h +++ b/include/vde/curves/bspline_surface.h @@ -14,6 +14,10 @@ public: int degree_u, int degree_v); [[nodiscard]] Point3D evaluate(double u, double v) const; + [[nodiscard]] Vector3D derivative_u(double u, double v) const; + [[nodiscard]] Vector3D derivative_v(double u, double v) const; + [[nodiscard]] Vector3D normal(double u, double v) const; + [[nodiscard]] int degree_u() const { return degree_u_; } [[nodiscard]] int degree_v() const { return degree_v_; } [[nodiscard]] const std::vector& knots_u() const { return knots_u_; } diff --git a/include/vde/curves/tessellation.h b/include/vde/curves/tessellation.h index 7050275..aca5dc2 100644 --- a/include/vde/curves/tessellation.h +++ b/include/vde/curves/tessellation.h @@ -2,12 +2,13 @@ #include "vde/core/point.h" #include #include +#include namespace vde::curves { using core::Point3D; using core::Vector3D; -/// Tessellate a parametric surface into a triangle mesh +/// Tessellate a parametric surface into a uniform triangle mesh /// @param eval Function(u,v) -> Point3D /// @param res_u, res_v Number of samples in u/v directions /// @return Pairs of (vertices, triangle indices) @@ -15,4 +16,20 @@ std::pair, std::vector>> tessellate(const std::function& eval, int res_u, int res_v); +/// Adaptive tessellation based on normal-angle deviation +/// Recursively subdivides the parameter domain when the angular deviation +/// between corner normals exceeds the threshold. +/// @param eval Function(u,v) -> Point3D +/// @param normal_fn Function(u,v) -> Vector3D (unit normal) +/// @param min_res Minimum subdivisions per direction (1 = single quad) +/// @param max_depth Maximum recursion depth (limits total subdivisions) +/// @param angle_threshold_deg Maximum normal-angle deviation in degrees +/// @return Pairs of (vertices, triangle indices) +std::pair, std::vector>> +adaptive_tessellate( + const std::function& eval, + const std::function& normal_fn, + int min_res = 2, int max_depth = 6, + double angle_threshold_deg = 5.0); + } // namespace vde::curves diff --git a/include/vde/foundation/io_obj.h b/include/vde/foundation/io_obj.h index 260499f..5409d35 100644 --- a/include/vde/foundation/io_obj.h +++ b/include/vde/foundation/io_obj.h @@ -7,8 +7,26 @@ namespace vde::foundation { struct ObjMeshData { std::vector vertices; + std::vector texcoords; std::vector normals; - std::vector> faces; // vertex indices (1-based in file, 0-based here) + + /// Vertex indices per face (0-based internally). + /// Always populated regardless of face format. + std::vector> faces; + + /// Texcoord indices per face (0-based). Empty if the file had no vt data. + /// Same length as faces; each sub-vector same length as corresponding face. + std::vector> face_texcoords; + + /// Normal indices per face (0-based). Empty if the file had no vn data. + std::vector> face_normals; + + /// Material name per face group (appears before the faces that use it). + /// The i-th entry gives the material active for the i-th face. + std::vector face_materials; + + /// Helpers + [[nodiscard]] bool has_texcoords() const { return !texcoords.empty(); } }; ObjMeshData read_obj(const std::string& filepath); diff --git a/include/vde/foundation/io_stl.h b/include/vde/foundation/io_stl.h index 2274085..8bbe61d 100644 --- a/include/vde/foundation/io_stl.h +++ b/include/vde/foundation/io_stl.h @@ -10,7 +10,13 @@ struct StlTriangle { Point3D v0, v1, v2; }; +/// Auto-detect format (binary or ASCII) and read. std::vector read_stl(const std::string& filepath); + +/// Always write binary STL. void write_stl(const std::string& filepath, const std::vector& tris); +/// Write ASCII STL (human-readable). +void write_stl_ascii(const std::string& filepath, const std::vector& tris); + } // namespace vde::foundation diff --git a/include/vde/mesh/mesh_quality.h b/include/vde/mesh/mesh_quality.h index 74cb893..37fe0f0 100644 --- a/include/vde/mesh/mesh_quality.h +++ b/include/vde/mesh/mesh_quality.h @@ -1,12 +1,16 @@ #pragma once #include "vde/mesh/halfedge_mesh.h" #include +#include +#include namespace vde::mesh { using core::Point2D; using core::Point3D; using core::Vector3D; +// ── Triangle quality ──────────────────────────────── + struct MeshQuality { double min_angle_deg = 0; double max_angle_deg = 0; @@ -17,4 +21,76 @@ struct MeshQuality { MeshQuality evaluate_mesh_quality(const HalfedgeMesh& mesh); +// ── Tetrahedron quality ───────────────────────────── + +struct TetQuality { + double min_jacobian = 0.0; // scaled Jacobian (normalised, 0=b−ad … 1=regular) + double max_jacobian = 0.0; + double avg_jacobian = 0.0; + double min_dihedral_deg = 0.0; // minimum dihedral angle (degrees) + double max_dihedral_deg = 0.0; // maximum dihedral angle + double max_aspect_ratio = 0.0; // circumsphere-radius / shortest-edge + size_t degenerate_tets = 0; + size_t total_tets = 0; +}; + +/// Evaluate tetrahedron quality from a triangle mesh. +/// This function interprets the mesh as a tetrahedral mesh where +/// each face belongs to a tet defined by the face vertices + face centroid +/// (suitable for closed triangle surfaces interpreted as volume boundaries). +TetQuality evaluate_tet_quality(const HalfedgeMesh& mesh); + +// ── Generic element quality ───────────────────────── + +enum class ElementType { Tri, Quad, Tet, Hex }; + +struct ElemQuality { + double min_jacobian = 0.0; + double max_jacobian = 0.0; + double avg_jacobian = 0.0; + size_t total_elements = 0; + size_t degenerate = 0; +}; + +/// Dispatch by element type; implemented for Tri and Tet (Quad/Hex return stub). +ElemQuality evaluate_element_quality(const HalfedgeMesh& mesh, ElementType type); + +// ── Quality report ────────────────────────────────── + +struct QualityReport { + std::string element_type; // "triangle" / "tetrahedron" + size_t total = 0; + size_t degenerate = 0; + + // Binned histogram (quality ∈ [0,1]) + static constexpr int kBins = 10; + std::array histogram{}; + + // Grade counts (A=excellent B=good C=acceptable D=poor F=fail) + size_t grade_A = 0, grade_B = 0, grade_C = 0, grade_D = 0, grade_F = 0; + + double min_quality = 0.0; + double max_quality = 0.0; + double avg_quality = 0.0; + double stddev_quality = 0.0; +}; + +/// Generate a full quality report from a triangle surface mesh. +QualityReport mesh_quality_report(const HalfedgeMesh& mesh); + +/// Generate a full quality report from a tetrahedral mesh. +/// The mesh is treated as having 4-vertex faces (tets) in the +/// same index order as the surface-mesh face loop. +QualityReport mesh_quality_report_tet(const HalfedgeMesh& mesh); + +// ── Utility ───────────────────────────────────────── + +/// Compute the scaled Jacobian (determinant / (product of edge-lengths)) +/// for a triangle; returns [−1,1], 1 = equilateral. +double tri_scaled_jacobian(const Point3D& a, const Point3D& b, const Point3D& c); + +/// Compute the scaled Jacobian for a tetrahedron; returns [0,1]. +double tet_scaled_jacobian(const Point3D& a, const Point3D& b, + const Point3D& c, const Point3D& d); + } // namespace vde::mesh diff --git a/include/vde/mesh/mesh_smooth.h b/include/vde/mesh/mesh_smooth.h index c251385..b15fba0 100644 --- a/include/vde/mesh/mesh_smooth.h +++ b/include/vde/mesh/mesh_smooth.h @@ -6,15 +6,38 @@ using core::Point2D; using core::Point3D; using core::Vector3D; -enum class SmoothMethod { Laplacian, Taubin }; +enum class SmoothMethod { Laplacian, Taubin, HCLaplacian, Bilateral }; struct SmoothOptions { int iterations = 10; - double lambda = 0.5; // Laplacian weight - double mu = -0.53; // Taubin shrink (mu < -lambda for Taubin) + double lambda = 0.5; // Laplacian weight + double mu = -0.53; // Taubin shrink (mu < -lambda for Taubin) + + // HC Laplacian parameters + double hc_alpha = 0.0; // push-back strength (0 = default auto) + double hc_beta = 0.5; // push-forward blend + + // Bilateral parameters + double bilateral_sigma_c = 0.0; // spatial sigma (0 = auto from avg edge length) + double bilateral_sigma_n = 0.3; // normal sigma (radians) + int bilateral_iters = 5; + SmoothMethod method = SmoothMethod::Laplacian; }; +/// Generic dispatcher — delegates to specific smooth functions HalfedgeMesh smooth_mesh(const HalfedgeMesh& mesh, const SmoothOptions& opts = {}); +/// Standard Laplacian smoothing +HalfedgeMesh smooth_laplacian(const HalfedgeMesh& mesh, int iterations, double lambda); + +/// Taubin λ|μ smoothing (volume-preserving) +HalfedgeMesh smooth_taubin(const HalfedgeMesh& mesh, int iterations, double lambda, double mu); + +/// HC Laplacian (Humphrey's Classes) — two-pass with push-back to preserve volume +HalfedgeMesh smooth_hc_laplacian(const HalfedgeMesh& mesh, const SmoothOptions& opts = {}); + +/// Bilateral mesh filtering — normal-weighted, preserves sharp features +HalfedgeMesh smooth_bilateral(const HalfedgeMesh& mesh, const SmoothOptions& opts = {}); + } // namespace vde::mesh diff --git a/include/vde/spatial/r_tree.h b/include/vde/spatial/r_tree.h index 0182bd3..01e5f11 100644 --- a/include/vde/spatial/r_tree.h +++ b/include/vde/spatial/r_tree.h @@ -1,6 +1,8 @@ #pragma once #include "vde/spatial/spatial_index.h" #include "vde/core/aabb.h" +#include +#include namespace vde::spatial { using core::Point3D; @@ -8,6 +10,15 @@ using core::Vector3D; using core::Ray3Dd; using core::Triangle3D; +template +struct RTreeNode { + AABB3D bbox; + bool is_leaf = true; + std::vector children; // internal: indices into nodes_ + std::vector item_indices; // leaf: indices into items_ + size_t parent = static_cast(-1); +}; + template class RTree : public SpatialIndex { public: @@ -19,9 +30,16 @@ public: std::vector query_ray(const Ray3Dd& ray) const override; void clear() override; size_t size() const override { return count_; } + size_t node_count() const { return nodes_.size(); } + private: + void query_range_recursive(size_t node_idx, const AABB3D& range, + std::vector& result) const; + std::vector items_; - std::vector bounds_; + std::vector item_bounds_; + std::vector> nodes_; + size_t root_idx_ = 0; size_t count_ = 0; }; diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8e94200..56cf214 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -20,8 +20,6 @@ target_link_libraries(vde_foundation # ── core ─────────────────────────────────────────── add_library(vde_core STATIC - core/point.cpp - core/transform.cpp core/polygon.cpp core/voronoi.cpp core/distance.cpp diff --git a/src/boolean/boolean_mesh.cpp b/src/boolean/boolean_mesh.cpp index e68fcfc..3252a38 100644 --- a/src/boolean/boolean_mesh.cpp +++ b/src/boolean/boolean_mesh.cpp @@ -1,6 +1,263 @@ #include "vde/boolean/boolean_mesh.h" +#include "vde/spatial/bvh.h" +#include "vde/core/triangle.h" +#include "vde/core/line.h" +#include +#include +#include +#include + namespace vde::boolean { -HalfedgeMesh mesh_boolean(const HalfedgeMesh&, const HalfedgeMesh&, BooleanOp) { +using core::Triangle3D; +using core::Ray3Dd; +using core::AABB3D; +using mesh::HalfedgeMesh; + +// ═══════════════════════════════════════════════════════════ +// Point-in-mesh classification via ray-casting +// ═══════════════════════════════════════════════════════════ + +static bool ray_tri_intersect( + const core::Point3D& origin, const core::Vector3D& dir, + const core::Point3D& v0, const core::Point3D& v1, const core::Point3D& v2) +{ + core::Vector3D e1 = v1 - v0; + core::Vector3D e2 = v2 - v0; + core::Vector3D h = dir.cross(e2); + double a = e1.dot(h); + if (std::abs(a) < 1e-12) return false; + + double f = 1.0 / a; + core::Vector3D s = origin - v0; + double u = f * s.dot(h); + if (u < 0.0 || u > 1.0) return false; + + core::Vector3D q = s.cross(e1); + double v = f * dir.dot(q); + if (v < 0.0 || u + v > 1.0) return false; + + double t = f * e2.dot(q); + return t > 1e-10; +} + +bool is_point_inside_mesh(const core::Point3D& p, const HalfedgeMesh& mesh) { + // Shoot a ray in a random-ish direction to avoid degenerate alignments + // Use a direction biased toward (1, 0.37, 0.53) to minimise edge/vertex grazing + core::Vector3D dir(1.0, 0.371, 0.539); + dir.normalize(); + + int hits = 0; + const size_t nf = mesh.num_faces(); + for (size_t fi = 0; fi < nf; ++fi) { + auto vis = mesh.face_vertices(static_cast(fi)); + if (vis.size() != 3) continue; + if (ray_tri_intersect(p, dir, + mesh.vertex(vis[0]), mesh.vertex(vis[1]), mesh.vertex(vis[2]))) + ++hits; + } + return (hits & 1) == 1; // odd → inside +} + +// ═══════════════════════════════════════════════════════════ +// BVH-accelerated point-in-mesh +// ═══════════════════════════════════════════════════════════ + +namespace { + +/// Build a BVH over `mesh` triangles. Returns an empty BVH if the mesh +/// has no triangles, in which case everything is "outside". +spatial::BVH build_bvh_for(const HalfedgeMesh& mesh) { + spatial::BVH bvh; + std::vector tris; + tris.reserve(mesh.num_faces()); + for (size_t fi = 0; fi < mesh.num_faces(); ++fi) { + auto vis = mesh.face_vertices(static_cast(fi)); + if (vis.size() == 3) + tris.emplace_back(mesh.vertex(vis[0]), + mesh.vertex(vis[1]), + mesh.vertex(vis[2])); + } + if (!tris.empty()) bvh.build(tris); + return bvh; +} + +/// Ray-cast using BVH; returns true if inside. +bool is_point_inside_bvh(const core::Point3D& p, + const spatial::BVH& bvh) { + if (bvh.size() == 0) return false; + + core::Vector3D dir(1.0, 0.371, 0.539); + dir.normalize(); + Ray3Dd ray(p, dir); + auto results = bvh.query_ray(ray); + int hits = 0; + for (const auto& tri : results) { + core::Vector3D e1 = tri.v(1) - tri.v(0); + core::Vector3D e2 = tri.v(2) - tri.v(0); + core::Vector3D h = dir.cross(e2); + double a = e1.dot(h); + if (std::abs(a) < 1e-12) continue; + double f = 1.0 / a; + core::Vector3D s = p - tri.v(0); + double u = f * s.dot(h); + if (u < 0.0 || u > 1.0) continue; + core::Vector3D q = s.cross(e1); + double v = f * dir.dot(q); + if (v < 0.0 || u + v > 1.0) continue; + double t = f * e2.dot(q); + if (t > 1e-10) ++hits; + } + return (hits & 1) == 1; +} + +} // anonymous namespace + +// ═══════════════════════════════════════════════════════════ +// Clip mesh — keep faces whose centroids are OUTSIDE b +// ═══════════════════════════════════════════════════════════ + +HalfedgeMesh clip_mesh(const HalfedgeMesh& a, const HalfedgeMesh& b, + bool invert_b) { + if (a.num_faces() == 0) return {}; + if (b.num_faces() == 0) { + // Empty boundary → everything is "outside" (inside inverted-empty) + return invert_b ? HalfedgeMesh{} : a; + } + + spatial::BVH bvh_b = build_bvh_for(b); + + std::vector out_verts; + std::vector> out_tris; + + for (size_t fi = 0; fi < a.num_faces(); ++fi) { + auto vis = a.face_vertices(static_cast(fi)); + if (vis.size() != 3) continue; + + core::Point3D v0 = a.vertex(vis[0]); + core::Point3D v1 = a.vertex(vis[1]); + core::Point3D v2 = a.vertex(vis[2]); + core::Point3D centroid = (v0 + v1 + v2) * (1.0 / 3.0); + + bool inside = is_point_inside_bvh(centroid, bvh_b); + if (invert_b) inside = !inside; + + // Keep the face if its centroid is OUTSIDE (b) + if (!inside) { + int idx = static_cast(out_verts.size()); + out_verts.push_back(v0); + out_verts.push_back(v1); + out_verts.push_back(v2); + // Preserve original winding (no flip) + out_tris.push_back({idx, idx + 1, idx + 2}); + } + } + + HalfedgeMesh result; + if (!out_verts.empty()) + result.build_from_triangles(out_verts, out_tris); + return result; +} + +// ═══════════════════════════════════════════════════════════ +// Mesh invert / merge helpers +// ═══════════════════════════════════════════════════════════ + +HalfedgeMesh invert_mesh(const HalfedgeMesh& m) { + // Flip face winding — reverse every triangle's vertex order + std::vector verts; + std::vector> tris; + verts.reserve(m.num_vertices()); + tris.reserve(m.num_faces()); + + for (size_t i = 0; i < m.num_vertices(); ++i) + verts.push_back(m.vertex(i)); + + for (size_t fi = 0; fi < m.num_faces(); ++fi) { + auto vis = m.face_vertices(static_cast(fi)); + if (vis.size() != 3) continue; + // Reverse winding: (v0, v1, v2) → (v0, v2, v1) + tris.push_back({vis[0], vis[2], vis[1]}); + } + + HalfedgeMesh result; + if (!verts.empty() && !tris.empty()) + result.build_from_triangles(verts, tris); + return result; +} + +HalfedgeMesh merge_meshes(const HalfedgeMesh& a, const HalfedgeMesh& b) { + if (a.num_vertices() == 0) return b; + if (b.num_vertices() == 0) return a; + + std::vector verts; + std::vector> tris; + + // Copy A + size_t offset_a = verts.size(); + for (size_t i = 0; i < a.num_vertices(); ++i) + verts.push_back(a.vertex(i)); + for (size_t fi = 0; fi < a.num_faces(); ++fi) { + auto vis = a.face_vertices(static_cast(fi)); + if (vis.size() != 3) continue; + tris.push_back({vis[0], vis[1], vis[2]}); + } + + // Copy B with vertex-index offset + int offset_b = static_cast(verts.size()); + for (size_t i = 0; i < b.num_vertices(); ++i) + verts.push_back(b.vertex(i)); + for (size_t fi = 0; fi < b.num_faces(); ++fi) { + auto vis = b.face_vertices(static_cast(fi)); + if (vis.size() != 3) continue; + tris.push_back({vis[0] + offset_b, + vis[1] + offset_b, + vis[2] + offset_b}); + } + + HalfedgeMesh result; + if (!verts.empty()) + result.build_from_triangles(verts, tris); + return result; +} + +// ═══════════════════════════════════════════════════════════ +// Public boolean operations +// ═══════════════════════════════════════════════════════════ + +HalfedgeMesh mesh_union(const HalfedgeMesh& a, const HalfedgeMesh& b) { + // Keep faces of each mesh that lie OUTSIDE the other + HalfedgeMesh a_out = clip_mesh(a, b, false); + HalfedgeMesh b_out = clip_mesh(b, a, false); + return merge_meshes(a_out, b_out); +} + +HalfedgeMesh mesh_intersection(const HalfedgeMesh& a, const HalfedgeMesh& b) { + // Keep faces whose centroids are INSIDE the other mesh. + // clip(…, invert=true) keeps faces OUTSIDE the inverted volume, + // i.e. INSIDE the original volume. + HalfedgeMesh a_in = clip_mesh(a, b, /*invert_b=*/true); + HalfedgeMesh b_in = clip_mesh(b, a, /*invert_b=*/true); + return merge_meshes(a_in, b_in); +} + +HalfedgeMesh mesh_difference(const HalfedgeMesh& a, const HalfedgeMesh& b) { + // A − B: keep A-faces that are OUTSIDE B + return clip_mesh(a, b, false); +} + +HalfedgeMesh mesh_sym_diff(const HalfedgeMesh& a, const HalfedgeMesh& b) { + return merge_meshes(mesh_difference(a, b), mesh_difference(b, a)); +} + +HalfedgeMesh mesh_boolean(const HalfedgeMesh& a, const HalfedgeMesh& b, + BooleanOp op) { + switch (op) { + case BooleanOp::Union: return mesh_union(a, b); + case BooleanOp::Intersection: return mesh_intersection(a, b); + case BooleanOp::Difference: return mesh_difference(a, b); + case BooleanOp::SymDiff: return mesh_sym_diff(a, b); + } return {}; } -} + +} // namespace vde::boolean diff --git a/src/collision/ray_intersect.cpp b/src/collision/ray_intersect.cpp index 6ed08f2..74e387e 100644 --- a/src/collision/ray_intersect.cpp +++ b/src/collision/ray_intersect.cpp @@ -1,5 +1,6 @@ #include "vde/collision/ray_intersect.h" #include +#include namespace vde::collision { using core::Point3D; @@ -7,6 +8,9 @@ using core::Vector3D; using core::Triangle3D; using core::Ray3Dd; using core::Segment3Dd; +using core::AABB3D; + +// ── Ray-Triangle (Möller-Trumbore) ─────────────────────────────────── std::optional ray_triangle_intersect( const Point3D& origin, const Vector3D& dir, const Triangle3D& tri) { @@ -28,4 +32,88 @@ std::optional ray_triangle_intersect( return RayTriResult{t, u, v, origin + dir * t}; } +// ── Ray-AABB (slab method) ─────────────────────────────────────────── + +bool ray_aabb_intersect(const Ray3Dd& ray, const AABB3D& box, + double& tmin_out, double& tmax_out) { + double tmin = -std::numeric_limits::infinity(); + double tmax = std::numeric_limits::infinity(); + + const Point3D& o = ray.origin(); + const Vector3D& d = ray.direction(); + const Point3D& bmin = box.min(); + const Point3D& bmax = box.max(); + + for (int i = 0; i < 3; ++i) { + if (std::abs(d[i]) < 1e-12) { + // Ray parallel to this slab pair + if (o[i] < bmin[i] || o[i] > bmax[i]) return false; + } else { + double inv_d = 1.0 / d[i]; + double t1 = (bmin[i] - o[i]) * inv_d; + double t2 = (bmax[i] - o[i]) * inv_d; + if (t1 > t2) std::swap(t1, t2); + tmin = std::max(tmin, t1); + tmax = std::min(tmax, t2); + if (tmin > tmax) return false; + } + } + + tmin_out = tmin; + tmax_out = tmax; + return tmax >= 0.0; +} + +// ── Ray-Sphere ─────────────────────────────────────────────────────── + +std::optional ray_sphere_intersect( + const Point3D& origin, const Vector3D& dir, + const Point3D& center, double radius) { + Vector3D oc = origin - center; + double a = dir.dot(dir); // = 1.0 for normalized dir + double b = 2.0 * oc.dot(dir); + double c = oc.dot(oc) - radius * radius; + double disc = b * b - 4.0 * a * c; + + if (disc < 0.0) return std::nullopt; + + double sqrt_disc = std::sqrt(disc); + double t = (-b - sqrt_disc) / (2.0 * a); + if (t < 0.0) + t = (-b + sqrt_disc) / (2.0 * a); + if (t < 0.0) return std::nullopt; + + return RaySphereResult{t, origin + dir * t}; +} + +// ── Ray-Plane ──────────────────────────────────────────────────────── + +std::optional ray_plane_intersect( + const Point3D& origin, const Vector3D& dir, + const Point3D& plane_point, const Vector3D& plane_normal) { + double denom = dir.dot(plane_normal); + if (std::abs(denom) < 1e-12) return std::nullopt; + double t = (plane_point - origin).dot(plane_normal) / denom; + if (t < 0.0) return std::nullopt; + return t; +} + +// ── Ray-Mesh (closest hit over all triangles) ──────────────────────── + +std::optional ray_mesh_intersect( + const Point3D& origin, const Vector3D& dir, + const std::vector& triangles) { + std::optional closest; + double closest_t = std::numeric_limits::infinity(); + + for (const auto& tri : triangles) { + auto hit = ray_triangle_intersect(origin, dir, tri); + if (hit && hit->t < closest_t) { + closest_t = hit->t; + closest = hit; + } + } + return closest; +} + } // namespace vde::collision diff --git a/src/collision/tri_intersect.cpp b/src/collision/tri_intersect.cpp index 58f51a5..ee2f630 100644 --- a/src/collision/tri_intersect.cpp +++ b/src/collision/tri_intersect.cpp @@ -1,25 +1,58 @@ #include "vde/collision/tri_intersect.h" #include +#include +#include +#include namespace vde::collision { using core::Point3D; using core::Vector3D; using core::Triangle3D; +using core::AABB3D; -// Separating axis test for two triangles -static bool on_opposite_sides(const Point3D& p, const Point3D& q, const Point3D& a, const Point3D& b) { - Vector3D pq = q - p, pa = a - p, pb = b - p; - return pa.dot(pq) * pb.dot(pq) < 0; +// ── helpers ─────────────────────────────────────────────────────────── + +namespace { + +/// Barycentric point-in-triangle test (3D, projection onto triangle plane). +bool point_in_tri(const Point3D& p, const Triangle3D& tri, double eps = 1e-10) { + Vector3D v0 = tri.v(2) - tri.v(0); + Vector3D v1 = tri.v(1) - tri.v(0); + Vector3D v2 = p - tri.v(0); + double d00 = v0.dot(v0), d01 = v0.dot(v1), d11 = v1.dot(v1); + double d20 = v2.dot(v0), d21 = v2.dot(v1); + double denom = d00 * d11 - d01 * d01; + if (std::abs(denom) < eps) return false; + double v = (d11 * d20 - d01 * d21) / denom; + double w = (d00 * d21 - d01 * d20) / denom; + return (v >= -eps) && (w >= -eps) && (v + w <= 1.0 + eps); } +/// Segment–plane intersection. Returns true and sets `out` if the +/// segment crosses the plane; result is clamped to [p0, p1]. +bool segment_plane_hit(const Point3D& p0, const Point3D& p1, + const Point3D& plane_pt, const Vector3D& plane_n, + Point3D& out) { + Vector3D dir = p1 - p0; + double denom = dir.dot(plane_n); + if (std::abs(denom) < 1e-12) return false; + double t = (plane_pt - p0).dot(plane_n) / denom; + if (t < -1e-12 || t > 1.0 + 1e-12) return false; + t = std::clamp(t, 0.0, 1.0); + out = p0 + dir * t; + return true; +} + +} // anonymous namespace + +// ── tri_tri_intersect (separating-axis, existing) ──────────────────── + bool tri_tri_intersect(const Triangle3D& t1, const Triangle3D& t2) { - // Compute plane of t2 Vector3D n2 = t2.normal(); double d1_0 = n2.dot(t1.v(0) - t2.v(0)); double d1_1 = n2.dot(t1.v(1) - t2.v(0)); double d1_2 = n2.dot(t1.v(2) - t2.v(0)); - // If all t1 vertices are on same side of t2 plane, no intersection if ((d1_0 > 1e-10 && d1_1 > 1e-10 && d1_2 > 1e-10) || (d1_0 < -1e-10 && d1_1 < -1e-10 && d1_2 < -1e-10)) return false; @@ -33,11 +66,9 @@ bool tri_tri_intersect(const Triangle3D& t1, const Triangle3D& t2) { (d2_0 < -1e-10 && d2_1 < -1e-10 && d2_2 < -1e-10)) return false; - // Compute intersection line L = n1 × n2 Vector3D L = n1.cross(n2); - if (L.norm() < 1e-12) return false; // Parallel planes + if (L.squaredNorm() < 1e-12) return false; - // Project vertices onto L and check interval overlap double t1_proj[3], t2_proj[3]; for (int i = 0; i < 3; ++i) t1_proj[i] = L.dot(t1.v(i)); for (int i = 0; i < 3; ++i) t2_proj[i] = L.dot(t2.v(i)); @@ -50,4 +81,139 @@ bool tri_tri_intersect(const Triangle3D& t1, const Triangle3D& t2) { return !(t1_max < t2_min || t2_max < t1_min); } +// ── tri_tri_intersect_detailed ─────────────────────────────────────── + +TriTriIntersection tri_tri_intersect_detailed(const Triangle3D& t1, + const Triangle3D& t2) { + Vector3D n1 = t1.normal(); + Vector3D n2 = t2.normal(); + Vector3D L = n1.cross(n2); + + if (L.squaredNorm() < 1e-12) { + // Coplanar — degenerate; quick check then bail with no segment + return {tri_tri_intersect(t1, t2), Point3D::Zero(), Point3D::Zero()}; + } + L.normalize(); + + // Collect all intersection points: edges of t1 ∩ plane(t2) + + // edges of t2 ∩ plane(t1), filtered to lie inside the other triangle. + std::vector pts; + + // t1 edges vs plane(t2) + for (int i = 0; i < 3; ++i) { + Point3D hit; + if (segment_plane_hit(t1.v(i), t1.v((i + 1) % 3), + t2.v(0), n2, hit)) { + if (point_in_tri(hit, t1) && point_in_tri(hit, t2)) + pts.push_back(hit); + } + } + + // t2 edges vs plane(t1) + for (int i = 0; i < 3; ++i) { + Point3D hit; + if (segment_plane_hit(t2.v(i), t2.v((i + 1) % 3), + t1.v(0), n1, hit)) { + if (point_in_tri(hit, t1) && point_in_tri(hit, t2)) + pts.push_back(hit); + } + } + + // Also include vertices that lie exactly on the other's plane + auto add_vertex_if_coplanar = [&](const Point3D& v, const Triangle3D& other, + const Vector3D& other_n) { + if (std::abs(other_n.dot(v - other.v(0))) < 1e-10) { + if (point_in_tri(v, other)) pts.push_back(v); + } + }; + for (int i = 0; i < 3; ++i) add_vertex_if_coplanar(t1.v(i), t2, n2); + for (int i = 0; i < 3; ++i) add_vertex_if_coplanar(t2.v(i), t1, n1); + + // Deduplicate by projection onto L, then sort + if (pts.size() < 2) { + return {!pts.empty(), pts.empty() ? Point3D::Zero() : pts[0], + pts.empty() ? Point3D::Zero() : pts[0]}; + } + + // Sort by projection onto L; collapse near-duplicates + std::sort(pts.begin(), pts.end(), + [&L](const Point3D& a, const Point3D& b) { + return L.dot(a) < L.dot(b); + }); + + std::vector uniq; + uniq.push_back(pts[0]); + for (size_t i = 1; i < pts.size(); ++i) { + if (L.dot(pts[i] - uniq.back()) > 1e-10) + uniq.push_back(pts[i]); + } + + if (uniq.size() < 2) { + return {true, uniq[0], uniq[0]}; + } + + return {true, uniq.front(), uniq.back()}; +} + +// ── tri_aabb_overlap (SAT with 13 separating axes) ─────────────────── + +bool tri_aabb_overlap(const Triangle3D& tri, const AABB3D& box) { + // Translate to box-center frame → box is symmetric about origin + Point3D center = box.center(); + Vector3D half = box.extent() * 0.5; + + Point3D v[3] = { + tri.v(0) - center, + tri.v(1) - center, + tri.v(2) - center + }; + + Vector3D e[3] = {v[1] - v[0], v[2] - v[1], v[0] - v[2]}; + Vector3D axes[3] = { + Vector3D(1, 0, 0), + Vector3D(0, 1, 0), + Vector3D(0, 0, 1) + }; + + // Build the 13 axes lazily: test as we go + auto test_axis = [&](const Vector3D& axis) -> bool { + double len2 = axis.squaredNorm(); + if (len2 < 1e-12) return true; // degenerate, skip + Vector3D d = axis / std::sqrt(len2); + + // Project triangle + double tri_min = d.dot(v[0]); + double tri_max = tri_min; + for (int i = 1; i < 3; ++i) { + double p = d.dot(v[i]); + tri_min = std::min(tri_min, p); + tri_max = std::max(tri_max, p); + } + + // Project box (half-extent onto axis) + double box_r = std::abs(d.x()) * half.x() + + std::abs(d.y()) * half.y() + + std::abs(d.z()) * half.z(); + + return !(tri_min > box_r || tri_max < -box_r); + }; + + // 1-3: AABB face normals + for (int i = 0; i < 3; ++i) + if (!test_axis(axes[i])) return false; + + // 4: triangle normal + Vector3D tri_n = e[0].cross(e[1]); + if (!test_axis(tri_n)) return false; + + // 5-13: cross products of triangle edges × AABB axes + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + if (!test_axis(e[i].cross(axes[j]))) return false; + } + } + + return true; +} + } // namespace vde::collision diff --git a/src/core/distance.cpp b/src/core/distance.cpp index 062736e..e5f6d97 100644 --- a/src/core/distance.cpp +++ b/src/core/distance.cpp @@ -1,8 +1,13 @@ #include "vde/core/distance.h" #include +#include namespace vde::core { +// ═══════════════════════════════════════════════ +// Point distances +// ═══════════════════════════════════════════════ + double distance(const Point3D& a, const Point3D& b) { return (a - b).norm(); } double distance(const Point3D& p, const Line3D& line) { @@ -11,11 +16,7 @@ double distance(const Point3D& p, const Line3D& line) { } double distance(const Point3D& p, const Segment3D& seg) { - Vector3D ab = seg.p1() - seg.p0(); - Vector3D ap = p - seg.p0(); - double t = ap.dot(ab) / ab.squaredNorm(); - t = std::clamp(t, 0.0, 1.0); - return (seg.p0() + t * ab - p).norm(); + return (closest_point(p, seg) - p).norm(); } double distance(const Point3D& p, const Plane& plane) { @@ -26,15 +27,143 @@ double distance(const Point3D& p, const Triangle& tri) { return (p - closest_point(p, tri)).norm(); } +double distance(const Point3D& p, const AABB& box) { + return (p - closest_point(p, box)).norm(); +} + +// ═══════════════════════════════════════════════ +// Feature-pair distances +// ═══════════════════════════════════════════════ + +double distance(const Segment3D& a, const Segment3D& b) { + // Closest points between two segments via the classic + // "clamp to [0,1] then clamp the other" approach. + Vector3D d1 = a.p1() - a.p0(); + Vector3D d2 = b.p1() - b.p0(); + Vector3D r = a.p0() - b.p0(); + + double a_len2 = d1.squaredNorm(); + double b_len2 = d2.squaredNorm(); + double ab = d1.dot(d2); + double ar = d1.dot(r); + double br = d2.dot(r); + + double denom = a_len2 * b_len2 - ab * ab; + double s, t; + + if (denom < 1e-12) { + // (Nearly) parallel segments + s = 0.0; + t = (ab > 0.0) ? br / b_len2 : 0.0; + t = std::clamp(t, 0.0, 1.0); + } else { + s = (ab * br - b_len2 * ar) / denom; + s = std::clamp(s, 0.0, 1.0); + t = (ab * s + br) / b_len2; + } + + // Clamp t then re-clamp s (standard "region" logic) + if (t < 0.0) { + t = 0.0; + s = std::clamp(-ar / a_len2, 0.0, 1.0); + } else if (t > 1.0) { + t = 1.0; + s = std::clamp((ab - ar) / a_len2, 0.0, 1.0); + } + + Point3D ca = a.p0() + s * d1; + Point3D cb = b.p0() + t * d2; + return (ca - cb).norm(); +} + +double distance(const Line3D& a, const Line3D& b) { + Vector3D u = a.direction(); + Vector3D v = b.direction(); + Vector3D w = a.origin() - b.origin(); + + double uv = u.dot(v); + double uu = u.squaredNorm(); // == 1 (normalized) + double vv = v.squaredNorm(); // == 1 + double wu = w.dot(u); + double wv = w.dot(v); + + double denom = uu * vv - uv * uv; + if (std::abs(denom) < 1e-12) { + // Parallel → distance from origin of a to line b + return (w - w.dot(v) * v).norm(); + } + double s = (uv * wv - vv * wu) / denom; + double t = (uu * wv - uv * wu) / denom; + + Point3D pa = a.origin() + s * u; + Point3D pb = b.origin() + t * v; + return (pa - pb).norm(); +} + +double distance(const AABB& a, const AABB& b) { + // Squared distance between two AABBs (0 when overlapping) + double sq = 0.0; + for (int i = 0; i < 3; ++i) { + if (a.max()(i) < b.min()(i)) { + double gap = b.min()(i) - a.max()(i); + sq += gap * gap; + } else if (b.max()(i) < a.min()(i)) { + double gap = a.min()(i) - b.max()(i); + sq += gap * gap; + } + // else: overlapping along this axis → contributes 0 + } + return std::sqrt(sq); +} + +double distance(const Triangle& a, const Triangle& b) { + // Minimum distance between two triangles is the minimum over: + // 6 vertex-face distances + 9 edge-edge distances + double dmin = std::numeric_limits::max(); + + // Vertex of a ↔ face of b + for (int i = 0; i < 3; ++i) { + dmin = std::min(dmin, distance(a.v(i), b)); + } + // Vertex of b ↔ face of a + for (int i = 0; i < 3; ++i) { + dmin = std::min(dmin, distance(b.v(i), a)); + } + // Edge of a ↔ edge of b + for (int i = 0; i < 3; ++i) { + Segment3D ea(a.v(i), a.v((i + 1) % 3)); + for (int j = 0; j < 3; ++j) { + Segment3D eb(b.v(j), b.v((j + 1) % 3)); + dmin = std::min(dmin, distance(ea, eb)); + } + } + return dmin; +} + +// ═══════════════════════════════════════════════ +// Closest-point queries +// ═══════════════════════════════════════════════ + Point3D closest_point(const Point3D& p, const Triangle& tri) { Vector3D v0 = tri.v(2) - tri.v(0), v1 = tri.v(1) - tri.v(0), v2 = p - tri.v(0); double d00 = v0.dot(v0), d01 = v0.dot(v1), d11 = v1.dot(v1); double d20 = v2.dot(v0), d21 = v2.dot(v1); double denom = d00 * d11 - d01 * d01; + if (std::abs(denom) < 1e-14) { + // Degenerate triangle – fall back to closest vertex + int best = 0; + double best_d2 = (p - tri.v(0)).squaredNorm(); + for (int i = 1; i < 3; ++i) { + double d2 = (p - tri.v(i)).squaredNorm(); + if (d2 < best_d2) { best_d2 = d2; best = i; } + } + return tri.v(best); + } double v = (d11 * d20 - d01 * d21) / denom; double w = (d00 * d21 - d01 * d20) / denom; if (v >= 0 && w >= 0 && v + w <= 1) return tri.v(0) + v * v0 + w * v1; + // Fallback to closest point on edges double d0 = distance(p, Segment3D(tri.v(0), tri.v(1))); double d1 = distance(p, Segment3D(tri.v(1), tri.v(2))); @@ -48,9 +177,24 @@ Point3D closest_point(const Point3D& p, const Triangle& tri) { Point3D closest_point(const Point3D& p, const Segment3D& seg) { Vector3D ab = seg.p1() - seg.p0(); Vector3D ap = p - seg.p0(); - double t = ap.dot(ab) / ab.squaredNorm(); + double len2 = ab.squaredNorm(); + if (len2 < 1e-14) return seg.p0(); + double t = ap.dot(ab) / len2; t = std::clamp(t, 0.0, 1.0); return seg.p0() + t * ab; } +Point3D closest_point(const Point3D& p, const AABB& box) { + // Clamp each coordinate to [min, max] + return Point3D( + std::clamp(p.x(), box.min().x(), box.max().x()), + std::clamp(p.y(), box.min().y(), box.max().y()), + std::clamp(p.z(), box.min().z(), box.max().z()) + ); +} + +Point3D closest_point(const Point3D& p, const Plane& plane) { + return plane.project(p); +} + } // namespace vde::core diff --git a/src/core/point.cpp b/src/core/point.cpp deleted file mode 100644 index 46bb60f..0000000 --- a/src/core/point.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#include "vde/core/point.h" -namespace vde::core {} diff --git a/src/core/polygon.cpp b/src/core/polygon.cpp index 8e3d07c..6459183 100644 --- a/src/core/polygon.cpp +++ b/src/core/polygon.cpp @@ -1,8 +1,73 @@ #include "vde/core/polygon.h" #include +#include namespace vde::core { +// ═══════════════════════════════════════════════ +// Perimeter +// ═══════════════════════════════════════════════ + +double Polygon2D::perimeter() const { + if (vertices_.size() < 2) return 0.0; + double p = 0.0; + for (size_t i = 0; i < vertices_.size(); ++i) { + const auto& a = vertices_[i]; + const auto& b = vertices_[(i + 1) % vertices_.size()]; + p += (a - b).norm(); + } + return p; +} + +// ═══════════════════════════════════════════════ +// Centroid (area-weighted) +// ═══════════════════════════════════════════════ + +Point2D Polygon2D::centroid() const { + if (vertices_.size() < 3) { + if (vertices_.empty()) return Point2D(0, 0); + // Point or line → average of vertices + Point2D sum(0, 0); + for (const auto& v : vertices_) sum += v; + return sum / static_cast(vertices_.size()); + } + + double A = signed_area(); + if (std::abs(A) < 1e-14) { + // Degenerate → average + Point2D sum(0, 0); + for (const auto& v : vertices_) sum += v; + return sum / static_cast(vertices_.size()); + } + + double cx = 0.0, cy = 0.0; + for (size_t i = 0; i < vertices_.size(); ++i) { + const auto& a = vertices_[i]; + const auto& b = vertices_[(i + 1) % vertices_.size()]; + double cross = a.x() * b.y() - b.x() * a.y(); + cx += (a.x() + b.x()) * cross; + cy += (a.y() + b.y()) * cross; + } + double factor = 1.0 / (6.0 * A); + return Point2D(cx * factor, cy * factor); +} + +// ═══════════════════════════════════════════════ +// Bounding box +// ═══════════════════════════════════════════════ + +AABB Polygon2D::bounding_box() const { + AABB aabb; + for (const auto& v : vertices_) { + aabb.expand(Point3D(v.x(), v.y(), 0.0)); + } + return aabb; +} + +// ═══════════════════════════════════════════════ +// Signed area (unchanged) +// ═══════════════════════════════════════════════ + double Polygon2D::signed_area() const { if (vertices_.size() < 3) return 0.0; double area = 0.0; @@ -14,8 +79,11 @@ double Polygon2D::signed_area() const { return area * 0.5; } +// ═══════════════════════════════════════════════ +// Point-in-polygon (unchanged) +// ═══════════════════════════════════════════════ + bool Polygon2D::contains(const Point2D& p) const { - // Ray casting algorithm bool inside = false; size_t n = vertices_.size(); for (size_t i = 0, j = n - 1; i < n; j = i++) { @@ -29,4 +97,193 @@ bool Polygon2D::contains(const Point2D& p) const { return inside; } +// ═══════════════════════════════════════════════ +// Douglas-Peucker simplification +// ═══════════════════════════════════════════════ + +namespace { + +double perpendicular_distance(const Point2D& p, const Point2D& a, const Point2D& b) { + Vector2D ab = b - a; + double len2 = ab.squaredNorm(); + if (len2 < 1e-14) return (p - a).norm(); + double t = (p - a).dot(ab) / len2; + t = std::clamp(t, 0.0, 1.0); + Point2D proj = a + t * ab; + return (p - proj).norm(); +} + +void douglas_peucker_recursive(const std::vector& pts, size_t start, size_t end, + double tolerance, std::vector& keep) { + if (end <= start + 1) return; + + double dmax = 0.0; + size_t idx = start; + for (size_t i = start + 1; i < end; ++i) { + double d = perpendicular_distance(pts[i], pts[start], pts[end]); + if (d > dmax) { + dmax = d; + idx = i; + } + } + + if (dmax > tolerance) { + keep[idx] = true; + douglas_peucker_recursive(pts, start, idx, tolerance, keep); + douglas_peucker_recursive(pts, idx, end, tolerance, keep); + } +} + +} // anonymous namespace + +Polygon2D Polygon2D::simplify(double tolerance) const { + if (vertices_.size() < 3) return *this; + + // For a closed polygon, we treat it as-is (first/last are not duplicated here). + // Apply DP on the chain 0..n-1, keeping endpoints. + std::vector keep(vertices_.size(), false); + keep.front() = true; + keep.back() = true; + + douglas_peucker_recursive(vertices_, 0, vertices_.size() - 1, tolerance, keep); + + std::vector result; + for (size_t i = 0; i < vertices_.size(); ++i) { + if (keep[i]) result.push_back(vertices_[i]); + } + return Polygon2D(std::move(result)); +} + +// ═══════════════════════════════════════════════ +// Ear-clipping triangulation +// ═══════════════════════════════════════════════ + +namespace { + +double cross_2d(const Point2D& a, const Point2D& b, const Point2D& c) { + return (b.x() - a.x()) * (c.y() - a.y()) - (b.y() - a.y()) * (c.x() - a.x()); +} + +bool point_in_triangle(const Point2D& p, const Point2D& a, + const Point2D& b, const Point2D& c) { + double s1 = cross_2d(a, b, p); + double s2 = cross_2d(b, c, p); + double s3 = cross_2d(c, a, p); + bool has_neg = (s1 < -1e-12) || (s2 < -1e-12) || (s3 < -1e-12); + bool has_pos = (s1 > 1e-12) || (s2 > 1e-12) || (s3 > 1e-12); + return !(has_neg && has_pos); +} + +bool is_ear(const std::vector& poly, int i, int j, int k, + const std::vector& remaining) { + // Check convexity (CCW polygon → cross should be > 0) + if (cross_2d(poly[i], poly[j], poly[k]) <= 0.0) return false; + + // Check no other vertex lies inside triangle (i,j,k) + for (int m : remaining) { + if (m == i || m == j || m == k) continue; + if (point_in_triangle(poly[m], poly[i], poly[j], poly[k])) return false; + } + return true; +} + +} // anonymous namespace + +std::vector> Polygon2D::triangulate() const { + std::vector> tris; + if (vertices_.size() < 3) return tris; + + // Work on a copy; use 0-based indices + std::vector poly = vertices_; + + // Ensure CCW + if (signed_area() < 0) { + std::reverse(poly.begin(), poly.end()); + } + + int n = static_cast(poly.size()); + std::vector prev(n), next(n), remaining; + for (int i = 0; i < n; ++i) { + prev[i] = (i - 1 + n) % n; + next[i] = (i + 1) % n; + remaining.push_back(i); + } + + int i = 0; + int fail_count = 0; + while (remaining.size() > 3 && fail_count < n * 2) { + int a = prev[i]; + int b = i; + int c = next[i]; + if (is_ear(poly, a, b, c, remaining)) { + tris.push_back({a, b, c}); + // Remove vertex b + next[a] = c; + prev[c] = a; + auto it = std::find(remaining.begin(), remaining.end(), b); + if (it != remaining.end()) remaining.erase(it); + i = a; // rewind slightly + fail_count = 0; + } else { + i = next[i]; + ++fail_count; + } + } + + // Final triangle + if (remaining.size() == 3) { + tris.push_back({remaining[0], remaining[1], remaining[2]}); + } + + return tris; +} + +// ═══════════════════════════════════════════════ +// Andrew's monotone chain convex hull +// ═══════════════════════════════════════════════ + +Polygon2D Polygon2D::convex_hull_2d(const std::vector& points) { + if (points.size() <= 1) return Polygon2D(points); + + auto sorted = points; + std::sort(sorted.begin(), sorted.end(), + [](const Point2D& a, const Point2D& b) { + return a.x() < b.x() || (a.x() == b.x() && a.y() < b.y()); + }); + + // Build lower hull + std::vector hull; + for (const auto& p : sorted) { + while (hull.size() >= 2) { + const auto& a = hull[hull.size() - 2]; + const auto& b = hull.back(); + if (cross_2d(a, b, p) <= 0.0) + hull.pop_back(); + else + break; + } + hull.push_back(p); + } + + // Build upper hull + size_t lower_size = hull.size(); + for (auto it = sorted.rbegin(); it != sorted.rend(); ++it) { + const auto& p = *it; + while (hull.size() > lower_size) { + const auto& a = hull[hull.size() - 2]; + const auto& b = hull.back(); + if (cross_2d(a, b, p) <= 0.0) + hull.pop_back(); + else + break; + } + hull.push_back(p); + } + + // Remove duplicate last point (== first point) + if (hull.size() > 1) hull.pop_back(); + + return Polygon2D(std::move(hull)); +} + } // namespace vde::core diff --git a/src/core/transform.cpp b/src/core/transform.cpp deleted file mode 100644 index 51d9c19..0000000 --- a/src/core/transform.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#include "vde/core/transform.h" -namespace vde::core {} diff --git a/src/curves/bezier_surface.cpp b/src/curves/bezier_surface.cpp index 9a32108..99d1911 100644 --- a/src/curves/bezier_surface.cpp +++ b/src/curves/bezier_surface.cpp @@ -1,4 +1,6 @@ #include "vde/curves/bezier_surface.h" +#include +#include namespace vde::curves { @@ -21,16 +23,54 @@ Point3D BezierSurface::evaluate(double u, double v) const { return de_casteljau(v, row_results); } -Vector3D BezierSurface::derivative_u(double, double) const { - return Vector3D::Zero(); // TODO +Vector3D BezierSurface::derivative_u(double u, double v) const { + int mu = degree_u(); + int mv = degree_v(); + if (mu < 1) return Vector3D::Zero(); + + // For each row (fixed v-index), compute derivative control points: + // mu * (P_{i+1,j} - P_{i,j}) for i = 0..mu-1 + // Evaluate in u, then combine in v + std::vector col; + col.reserve(mv + 1); + for (const auto& row : cp_) { + std::vector dcp_u; + dcp_u.reserve(mu); + for (int i = 0; i < mu; ++i) + dcp_u.push_back(mu * (row[i + 1] - row[i])); + col.push_back(de_casteljau(u, dcp_u)); + } + return de_casteljau(v, col) - Point3D::Zero(); } -Vector3D BezierSurface::derivative_v(double, double) const { - return Vector3D::Zero(); // TODO +Vector3D BezierSurface::derivative_v(double u, double v) const { + int mu = degree_u(); + int mv = degree_v(); + if (mv < 1) return Vector3D::Zero(); + + // For each column (fixed u-index), compute derivative control points: + // mv * (P_{i,j+1} - P_{i,j}) for j = 0..mv-1 + // Evaluate in v, then combine in u + std::vector v_evals; + v_evals.reserve(mu + 1); + for (int i = 0; i <= mu; ++i) { + std::vector dcp_v; + dcp_v.reserve(mv); + for (int j = 0; j < mv; ++j) + dcp_v.push_back(mv * (cp_[j + 1][i] - cp_[j][i])); + v_evals.push_back(de_casteljau(v, dcp_v)); + } + return de_casteljau(u, v_evals) - Point3D::Zero(); } Vector3D BezierSurface::normal(double u, double v) const { - return derivative_u(u, v).cross(derivative_v(u, v)).normalized(); + Vector3D du = derivative_u(u, v); + Vector3D dv = derivative_v(u, v); + Vector3D n = du.cross(dv); + double nrm = n.norm(); + if (nrm > 1e-12) + return n / nrm; + return Vector3D::UnitZ(); } } // namespace vde::curves diff --git a/src/curves/bspline_curve.cpp b/src/curves/bspline_curve.cpp index ddde668..096f2c9 100644 --- a/src/curves/bspline_curve.cpp +++ b/src/curves/bspline_curve.cpp @@ -1,5 +1,6 @@ #include "vde/curves/bspline_curve.h" #include +#include namespace vde::curves { @@ -42,6 +43,67 @@ std::vector BSplineCurve::basis_functions(double t, int span) const { return N; } +std::vector BSplineCurve::basis_derivatives(double t, int span) const { + if (span < 0) span = find_span(t); + + // For degree 0, derivatives are zero + if (degree_ <= 0) + return std::vector(1, 0.0); + + // Build the ndu matrix: ndu[k][i] = N_{span-degree+k, i}(t) + // This is a full triangular storage of all intermediate basis values. + // ndu[j][r] stores the denominator for the recurrence when building column j, row r, + // and ndu[r][j] stores the triangular value N_{span-j+r, j}. + std::vector> ndu(degree_ + 1, + std::vector(degree_ + 1, 0.0)); + ndu[0][0] = 1.0; + + std::vector left(degree_ + 1), right(degree_ + 1); + + for (int j = 1; j <= degree_; ++j) { + left[j] = t - knots_[span + 1 - j]; + right[j] = knots_[span + j] - t; + double saved = 0.0; + for (int r = 0; r < j; ++r) { + ndu[j][r] = right[r + 1] + left[j - r]; // denominator + double temp = ndu[r][j - 1] / ndu[j][r]; + ndu[r][j] = saved + right[r + 1] * temp; + saved = left[j - r] * temp; + } + ndu[j][j] = saved; + } + + // Compute first derivatives using the formula: + // N'_{i,p} = p * ( N_{i,p-1} / (u_{i+p} - u_i) + // - N_{i+1,p-1} / (u_{i+p+1} - u_{i+1}) ) + // For i = span - degree + r (r = 0..degree), this becomes: + // ders[r] = degree * ( ndu[r-1][degree-1] / (u_{span+r} - u_{span-degree+r}) + // - ndu[r][degree-1] / (u_{span+r+1} - u_{span-degree+r+1}) ) + // with ndu[-1][*] = 0 and ndu[degree][degree-1] = 0. + + std::vector ders(degree_ + 1, 0.0); + + for (int r = 0; r <= degree_; ++r) { + double term1 = 0.0, term2 = 0.0; + + if (r > 0) { + double denom = knots_[span + r] - knots_[span - degree_ + r]; + if (std::abs(denom) > 1e-15) + term1 = ndu[r - 1][degree_ - 1] / denom; + } + + if (r < degree_) { + double denom = knots_[span + r + 1] - knots_[span - degree_ + r + 1]; + if (std::abs(denom) > 1e-15) + term2 = ndu[r][degree_ - 1] / denom; + } + + ders[r] = degree_ * (term1 - term2); + } + + return ders; +} + Point3D BSplineCurve::evaluate(double t) const { int span = find_span(t); auto N = basis_functions(t, span); diff --git a/src/curves/bspline_surface.cpp b/src/curves/bspline_surface.cpp index 38a5aaf..200ae52 100644 --- a/src/curves/bspline_surface.cpp +++ b/src/curves/bspline_surface.cpp @@ -1,4 +1,5 @@ #include "vde/curves/bspline_surface.h" +#include namespace vde::curves { @@ -23,4 +24,46 @@ Point3D BSplineSurface::evaluate(double u, double v) const { return result; } +Vector3D BSplineSurface::derivative_u(double u, double v) const { + BSplineCurve bs_u({}, knots_u_, degree_u_); + int span_u = bs_u.find_span(u); + auto dNu = bs_u.basis_derivatives(u, span_u); + + BSplineCurve bs_v({}, knots_v_, degree_v_); + int span_v = bs_v.find_span(v); + auto Nv = bs_v.basis_functions(v, span_v); + + Point3D result = Point3D::Zero(); + for (int i = 0; i <= degree_u_; ++i) + for (int j = 0; j <= degree_v_; ++j) + result += dNu[i] * Nv[j] * cp_[span_u - degree_u_ + i][span_v - degree_v_ + j]; + return result - Point3D::Zero(); +} + +Vector3D BSplineSurface::derivative_v(double u, double v) const { + BSplineCurve bs_u({}, knots_u_, degree_u_); + int span_u = bs_u.find_span(u); + auto Nu = bs_u.basis_functions(u, span_u); + + BSplineCurve bs_v({}, knots_v_, degree_v_); + int span_v = bs_v.find_span(v); + auto dNv = bs_v.basis_derivatives(v, span_v); + + Point3D result = Point3D::Zero(); + for (int i = 0; i <= degree_u_; ++i) + for (int j = 0; j <= degree_v_; ++j) + result += Nu[i] * dNv[j] * cp_[span_u - degree_u_ + i][span_v - degree_v_ + j]; + return result - Point3D::Zero(); +} + +Vector3D BSplineSurface::normal(double u, double v) const { + Vector3D du = derivative_u(u, v); + Vector3D dv = derivative_v(u, v); + Vector3D n = du.cross(dv); + double nrm = n.norm(); + if (nrm > 1e-12) + return n / nrm; + return Vector3D::UnitZ(); +} + } // namespace vde::curves diff --git a/src/curves/tessellation.cpp b/src/curves/tessellation.cpp index ef10833..0483355 100644 --- a/src/curves/tessellation.cpp +++ b/src/curves/tessellation.cpp @@ -1,7 +1,12 @@ #include "vde/curves/tessellation.h" +#include +#include +#include namespace vde::curves { +// --- Uniform tessellation --- + std::pair, std::vector>> tessellate(const std::function& eval, int res_u, int res_v) { @@ -31,4 +36,114 @@ tessellate(const std::function& eval, return {verts, tris}; } +// --- Adaptive tessellation --- + +namespace { + +/// Hash a (u,v) parameter pair for the vertex cache +struct ParamKey { + double u, v; + bool operator==(const ParamKey& o) const { + return std::abs(u - o.u) < 1e-12 && std::abs(v - o.v) < 1e-12; + } +}; + +struct ParamKeyHash { + std::size_t operator()(const ParamKey& k) const { + auto hu = std::hash{}(k.u); + auto hv = std::hash{}(k.v); + return hu ^ (hv + 0x9e3779b9 + (hu << 6) + (hu >> 2)); + } +}; + +/// Compute the angle (in degrees) between two vectors +inline double angle_deg(const Vector3D& a, const Vector3D& b) { + double dot = std::clamp(a.normalized().dot(b.normalized()), -1.0, 1.0); + return std::acos(dot) * 180.0 / M_PI; +} + +/// Recursive adaptive subdivision of a patch [u0,u1]×[v0,v1] +void subdivide( + double u0, double u1, double v0, double v1, + const std::function& eval, + const std::function& normal_fn, + int depth, int max_depth, double cos_threshold, + std::vector& verts, + std::vector>& tris, + std::unordered_map& vcache) +{ + // Evaluate normals at the four corners + Vector3D n00 = normal_fn(u0, v0); + Vector3D n10 = normal_fn(u1, v0); + Vector3D n01 = normal_fn(u0, v1); + Vector3D n11 = normal_fn(u1, v1); + + // Check maximum angular deviation + double max_angle = std::max({ + angle_deg(n00, n10), angle_deg(n00, n01), + angle_deg(n00, n11), angle_deg(n10, n01), + angle_deg(n10, n11), angle_deg(n01, n11) + }); + + // Stop subdividing if within tolerance or max depth reached + if (depth >= max_depth || std::cos(max_angle * M_PI / 180.0) >= cos_threshold) { + // Emit two triangles for this patch + auto get_idx = [&](double u, double v) -> int { + ParamKey k{u, v}; + auto it = vcache.find(k); + if (it != vcache.end()) return it->second; + int idx = static_cast(verts.size()); + verts.push_back(eval(u, v)); + vcache[k] = idx; + return idx; + }; + + int a = get_idx(u0, v0); + int b = get_idx(u1, v0); + int c = get_idx(u0, v1); + int d = get_idx(u1, v1); + tris.push_back({a, b, d}); + tris.push_back({a, d, c}); + return; + } + + double um = 0.5 * (u0 + u1); + double vm = 0.5 * (v0 + v1); + int next_depth = depth + 1; + + subdivide(u0, um, v0, vm, eval, normal_fn, next_depth, max_depth, cos_threshold, verts, tris, vcache); + subdivide(um, u1, v0, vm, eval, normal_fn, next_depth, max_depth, cos_threshold, verts, tris, vcache); + subdivide(u0, um, vm, v1, eval, normal_fn, next_depth, max_depth, cos_threshold, verts, tris, vcache); + subdivide(um, u1, vm, v1, eval, normal_fn, next_depth, max_depth, cos_threshold, verts, tris, vcache); +} + +} // anonymous namespace + +std::pair, std::vector>> +adaptive_tessellate( + const std::function& eval, + const std::function& normal_fn, + int min_res, int max_depth, double angle_threshold_deg) { + + std::vector verts; + std::vector> tris; + std::unordered_map vcache; + + double cos_threshold = std::cos(angle_threshold_deg * M_PI / 180.0); + + // Start with a uniform base grid at min_res, then subdivide each cell adaptively + for (int j = 0; j < min_res; ++j) { + double v0 = static_cast(j) / min_res; + double v1 = static_cast(j + 1) / min_res; + for (int i = 0; i < min_res; ++i) { + double u0 = static_cast(i) / min_res; + double u1 = static_cast(i + 1) / min_res; + subdivide(u0, u1, v0, v1, eval, normal_fn, 0, max_depth, + cos_threshold, verts, tris, vcache); + } + } + + return {verts, tris}; +} + } // namespace vde::curves diff --git a/src/foundation/io_obj.cpp b/src/foundation/io_obj.cpp index 52107ca..a41a6d9 100644 --- a/src/foundation/io_obj.cpp +++ b/src/foundation/io_obj.cpp @@ -10,29 +10,76 @@ ObjMeshData read_obj(const std::string& filepath) { if (!file) return data; std::string line; + std::string current_material; + while (std::getline(file, line)) { + // Trim trailing \r + if (!line.empty() && line.back() == '\r') line.pop_back(); + if (line.empty() || line[0] == '#') continue; std::istringstream iss(line); std::string token; iss >> token; + if (token == "v") { double x, y, z; iss >> x >> y >> z; data.vertices.emplace_back(x, y, z); - } else if (token == "vn") { + } + else if (token == "vt") { + double u, v, w = 0.0; + iss >> u >> v; + if (iss >> w) { /* 1D/3D texcoord — keep only uv */ } + data.texcoords.emplace_back(u, v); + } + else if (token == "vn") { double nx, ny, nz; iss >> nx >> ny >> nz; data.normals.emplace_back(nx, ny, nz); - } else if (token == "f") { - std::vector face; + } + else if (token == "usemtl") { + iss >> current_material; + } + else if (token == "f") { + std::vector face_v; + std::vector face_vt; + std::vector face_vn; + std::string vertex; while (iss >> vertex) { - std::istringstream viss(vertex); - std::string vi_str; - std::getline(viss, vi_str, '/'); - face.push_back(std::stoi(vi_str) - 1); + // Parse v/vt/vn, v//vn, v/vt formats + int v_idx = -1; + int vt_idx = -1; + int vn_idx = -1; + + // Count slashes + size_t slash1 = vertex.find('/'); + if (slash1 == std::string::npos) { + // "v" only + v_idx = std::stoi(vertex) - 1; + } else { + v_idx = std::stoi(vertex.substr(0, slash1)) - 1; + size_t slash2 = vertex.find('/', slash1 + 1); + if (slash2 == std::string::npos) { + // "v/vt" + vt_idx = std::stoi(vertex.substr(slash1 + 1)) - 1; + } else { + // "v/vt/vn" or "v//vn" + std::string mid = vertex.substr(slash1 + 1, slash2 - slash1 - 1); + if (!mid.empty()) vt_idx = std::stoi(mid) - 1; + vn_idx = std::stoi(vertex.substr(slash2 + 1)) - 1; + } + } + + face_v.push_back(v_idx); + face_vt.push_back(vt_idx); + face_vn.push_back(vn_idx); } - data.faces.push_back(std::move(face)); + + data.faces.push_back(std::move(face_v)); + data.face_texcoords.push_back(std::move(face_vt)); + data.face_normals.push_back(std::move(face_vn)); + data.face_materials.push_back(current_material); } } return data; @@ -40,13 +87,57 @@ ObjMeshData read_obj(const std::string& filepath) { void write_obj(const std::string& filepath, const ObjMeshData& data) { std::ofstream file(filepath); + if (!file) return; + + file.precision(12); + + // Vertices for (const auto& v : data.vertices) file << "v " << v.x() << " " << v.y() << " " << v.z() << "\n"; + + // Texture coordinates + for (const auto& vt : data.texcoords) + file << "vt " << vt.x() << " " << vt.y() << "\n"; + + // Normals + bool has_vt = data.has_texcoords(); + bool has_vn = !data.normals.empty(); + for (const auto& n : data.normals) file << "vn " << n.x() << " " << n.y() << " " << n.z() << "\n"; - for (const auto& f : data.faces) { + + // Faces with material tracking + std::string current_mtl; + for (size_t fi = 0; fi < data.faces.size(); ++fi) { + // Emit material change + if (fi < data.face_materials.size() && data.face_materials[fi] != current_mtl) { + current_mtl = data.face_materials[fi]; + if (!current_mtl.empty()) + file << "usemtl " << current_mtl << "\n"; + } + file << "f"; - for (int vi : f) file << " " << (vi + 1); + const auto& fv = data.faces[fi]; + for (size_t vi = 0; vi < fv.size(); ++vi) { + file << " " << (fv[vi] + 1); // 0-based → 1-based + + bool wrote_slash = false; + if (has_vt && fi < data.face_texcoords.size() && vi < data.face_texcoords[fi].size()) { + int vt = data.face_texcoords[fi][vi]; + if (vt >= 0) { + if (!wrote_slash) { file << "/"; wrote_slash = true; } + file << (vt + 1); + } + } + + if (has_vn && fi < data.face_normals.size() && vi < data.face_normals[fi].size()) { + int vn = data.face_normals[fi][vi]; + if (vn >= 0) { + if (!wrote_slash) file << "/"; + file << "/" << (vn + 1); + } + } + } file << "\n"; } } diff --git a/src/foundation/io_stl.cpp b/src/foundation/io_stl.cpp index 54fe738..cb29146 100644 --- a/src/foundation/io_stl.cpp +++ b/src/foundation/io_stl.cpp @@ -1,21 +1,49 @@ #include "vde/foundation/io_stl.h" #include -#include +#include +#include namespace vde::foundation { -std::vector read_stl(const std::string& filepath) { +// ── helpers ─────────────────────────────────────── + +namespace { + +bool looks_like_ascii_stl(const std::string& filepath) { + std::ifstream file(filepath); + if (!file) return false; + std::string line; + std::getline(file, line); + // Trim leading whitespace on the first token + std::istringstream iss(line); + std::string token; + iss >> token; + // Lowercase for case-insensitive check + for (auto& c : token) c = static_cast(std::tolower(c)); + return token == "solid"; +} + +std::vector read_stl_binary(const std::string& filepath) { std::vector tris; std::ifstream file(filepath, std::ios::binary); if (!file) return tris; - // Read binary STL char header[80]; file.read(header, 80); uint32_t count; file.read(reinterpret_cast(&count), 4); - tris.reserve(count); + // Sanity check: file size should match + file.seekg(0, std::ios::end); + auto fsize = file.tellg(); + file.seekg(84, std::ios::beg); + auto expected = static_cast(84 + 50 * count); + if (fsize < expected) { + // Not a valid binary STL + return {}; + } + + tris.reserve(count); for (uint32_t i = 0; i < count; ++i) { StlTriangle tri{}; file.read(reinterpret_cast(&tri.normal), 12); @@ -29,6 +57,64 @@ std::vector read_stl(const std::string& filepath) { return tris; } +std::vector read_stl_ascii(const std::string& filepath) { + std::vector tris; + std::ifstream file(filepath); + if (!file) return tris; + + std::string line; + StlTriangle current{}; + bool in_facet = false; + int vertex_count = 0; + + while (std::getline(file, line)) { + if (!line.empty() && line.back() == '\r') line.pop_back(); + + std::istringstream iss(line); + std::string token; + iss >> token; + + // Lowercase for case-insensitive matching + for (auto& c : token) c = static_cast(std::tolower(c)); + + if (token == "facet" || token == "facet") { + in_facet = true; + vertex_count = 0; + current = StlTriangle{}; + // Read "normal nx ny nz" + std::string normal_kw; + iss >> normal_kw; + iss >> current.normal.x() >> current.normal.y() >> current.normal.z(); + } + else if (token == "vertex" && in_facet) { + double x, y, z; + iss >> x >> y >> z; + if (vertex_count == 0) current.v0 = Point3D(x, y, z); + else if (vertex_count == 1) current.v1 = Point3D(x, y, z); + else if (vertex_count == 2) current.v2 = Point3D(x, y, z); + ++vertex_count; + } + else if (token == "endfacet") { + in_facet = false; + tris.push_back(current); + } + } + return tris; +} + +} // anonymous namespace + +// ── public API ──────────────────────────────────── + +std::vector read_stl(const std::string& filepath) { + if (looks_like_ascii_stl(filepath)) { + auto tris = read_stl_ascii(filepath); + if (!tris.empty()) return tris; + // Fall through – maybe "solid" is just a coincidental header in a binary STL + } + return read_stl_binary(filepath); +} + void write_stl(const std::string& filepath, const std::vector& tris) { std::ofstream file(filepath, std::ios::binary); char header[80] = {}; @@ -45,4 +131,21 @@ void write_stl(const std::string& filepath, const std::vector& tris } } +void write_stl_ascii(const std::string& filepath, const std::vector& tris) { + std::ofstream file(filepath); + file.precision(12); + file << "solid vde\n"; + for (const auto& tri : tris) { + file << " facet normal " + << tri.normal.x() << " " << tri.normal.y() << " " << tri.normal.z() << "\n"; + file << " outer loop\n"; + file << " vertex " << tri.v0.x() << " " << tri.v0.y() << " " << tri.v0.z() << "\n"; + file << " vertex " << tri.v1.x() << " " << tri.v1.y() << " " << tri.v1.z() << "\n"; + file << " vertex " << tri.v2.x() << " " << tri.v2.y() << " " << tri.v2.z() << "\n"; + file << " endloop\n"; + file << " endfacet\n"; + } + file << "endsolid vde\n"; +} + } // namespace vde::foundation diff --git a/src/mesh/mesh_quality.cpp b/src/mesh/mesh_quality.cpp index a00c461..82816ba 100644 --- a/src/mesh/mesh_quality.cpp +++ b/src/mesh/mesh_quality.cpp @@ -1,9 +1,16 @@ #include "vde/mesh/mesh_quality.h" +#include "vde/core/plane.h" #include #include +#include +#include namespace vde::mesh { +// ═══════════════════════════════════════════════════════════ +// Triangle quality (existing, kept for reference) +// ═══════════════════════════════════════════════════════════ + MeshQuality evaluate_mesh_quality(const HalfedgeMesh& mesh) { MeshQuality q{}; if (mesh.num_faces() == 0) return q; @@ -50,4 +57,393 @@ MeshQuality evaluate_mesh_quality(const HalfedgeMesh& mesh) { return q; } +// ═══════════════════════════════════════════════════════════ +// Utilities: scaled Jacobian for triangle / tetrahedron +// ═══════════════════════════════════════════════════════════ + +double tri_scaled_jacobian(const Point3D& a, const Point3D& b, const Point3D& c) { + Vector3D e1 = b - a; + Vector3D e2 = c - a; + double L1 = e1.norm(); + double L2 = e2.norm(); + double L3 = (c - b).norm(); + if (L1 < 1e-12 || L2 < 1e-12 || L3 < 1e-12) return 0.0; + // Area of equilateral triangle with same average edge-length + double avg_L = (L1 + L2 + L3) / 3.0; + double area_ref = (avg_L * avg_L) * std::sqrt(3.0) / 4.0; + double area = e1.cross(e2).norm() * 0.5; + return std::min(1.0, area / area_ref); +} + +double tet_scaled_jacobian(const Point3D& a, const Point3D& b, + const Point3D& c, const Point3D& d) { + // Edge vectors from vertex a + Vector3D e1 = b - a; + Vector3D e2 = c - a; + Vector3D e3 = d - a; + + // Determinant (6× signed volume) + double detJ = std::abs(e1.dot(e2.cross(e3))); + + // Product of edge-lengths (L2-norm product) + double prod = e1.norm() * e2.norm() * e3.norm(); + if (prod < 1e-12) return 0.0; + + // Normalise: regular tet has |detJ| / (edge-product) = √2/2 ≈ 0.707... + // Scale to [0,1] + return std::min(1.0, detJ / prod / 0.7071067811865476); +} + +// ═══════════════════════════════════════════════════════════ +// Tetrahedron quality evaluation +// ═══════════════════════════════════════════════════════════ + +static double tet_dihedral_angle(const Point3D& a, const Point3D& b, + const Point3D& c, const Point3D& d) { + // Compute all 6 dihedral angles; return the range [min, max] + // Face normals for each of the 4 faces + auto face_norm = [](const Point3D& p0, const Point3D& p1, const Point3D& p2) -> Vector3D { + return (p1 - p0).cross(p2 - p0).normalized(); + }; + + Vector3D n_abc = face_norm(a, b, c); // opposite d + Vector3D n_abd = face_norm(a, b, d); // opposite c + Vector3D n_acd = face_norm(a, c, d); // opposite b + Vector3D n_bcd = face_norm(b, c, d); // opposite a + + // Flip normals to point outward from centroid + Point3D centroid = (a + b + c + d) * 0.25; + if (n_abc.dot(a - centroid) < 0) n_abc = -n_abc; + if (n_abd.dot(a - centroid) < 0) n_abd = -n_abd; + if (n_acd.dot(a - centroid) < 0) n_acd = -n_acd; + if (n_bcd.dot(b - centroid) < 0) n_bcd = -n_bcd; + + auto dihedral = [](const Vector3D& n1, const Vector3D& n2) -> double { + double dot = std::clamp(n1.dot(n2), -1.0, 1.0); + return std::acos(dot) * 180.0 / M_PI; + }; + + // Return the min dihedral for this tet (worst angle) + double angle_ab = dihedral(n_abc, n_abd); // edge ab + double angle_ac = dihedral(n_abc, -n_acd); // edge ac + double angle_ad = dihedral(n_abd, -n_acd); // edge ad + double angle_bc = dihedral(n_abc, -n_bcd); // edge bc + double angle_bd = dihedral(n_abd, -n_bcd); // edge bd + double angle_cd = dihedral(n_acd, -n_bcd); // edge cd + + return std::min({angle_ab, angle_ac, angle_ad, angle_bc, angle_bd, angle_cd}); +} + +static double tet_aspect_ratio(const Point3D& a, const Point3D& b, + const Point3D& c, const Point3D& d) { + // Aspect ratio = circumsphere radius / shortest edge + double l_ab = (b - a).norm(); + double l_ac = (c - a).norm(); + double l_ad = (d - a).norm(); + double l_bc = (c - b).norm(); + double l_bd = (d - b).norm(); + double l_cd = (d - c).norm(); + double min_edge = std::min({l_ab, l_ac, l_ad, l_bc, l_bd, l_cd}); + if (min_edge < 1e-12) return 1e6; + + // Circumsphere radius via Cayley-Menger determinant + // R = √(D_circ / (2 * D_vol)) + auto sq = [](double x) { return x * x; }; + + double d_ab = sq(l_ab), d_ac = sq(l_ac), d_ad = sq(l_ad); + double d_bc = sq(l_bc), d_bd = sq(l_bd), d_cd = sq(l_cd); + + // Cayley-Menger determinant for circumradius + double M[5][5] = { + {0, 1, 1, 1, 1}, + {1, 0, d_ab, d_ac, d_ad}, + {1, d_ab, 0, d_bc, d_bd}, + {1, d_ac, d_bc, 0, d_cd}, + {1, d_ad, d_bd, d_cd, 0} + }; + + // 4×4 co-factor for volume: skip first row/col + double M_vol[4][4] = { + {0, d_ab, d_ac, d_ad}, + {d_ab, 0, d_bc, d_bd}, + {d_ac, d_bc, 0, d_cd}, + {d_ad, d_bd, d_cd, 0 } + }; + + // Compute determinant via simple 4×4 expansion + auto det4 = [](const double m[4][4]) -> double { + double a00 = m[0][0], a01 = m[0][1], a02 = m[0][2], a03 = m[0][3]; + double a10 = m[1][0], a11 = m[1][1], a12 = m[1][2], a13 = m[1][3]; + double a20 = m[2][0], a21 = m[2][1], a22 = m[2][2], a23 = m[2][3]; + double a30 = m[3][0], a31 = m[3][1], a32 = m[3][2], a33 = m[3][3]; + + return a00 * a11 * a22 * a33 + a00 * a12 * a23 * a31 + a00 * a13 * a21 * a32 + + a01 * a10 * a23 * a32 + a01 * a12 * a20 * a33 + a01 * a13 * a22 * a30 + + a02 * a10 * a21 * a33 + a02 * a11 * a23 * a30 + a02 * a13 * a20 * a31 + + a03 * a10 * a22 * a31 + a03 * a11 * a20 * a32 + a03 * a12 * a21 * a30 + - a00 * a11 * a23 * a32 - a00 * a12 * a21 * a33 - a00 * a13 * a22 * a31 + - a01 * a10 * a22 * a33 - a01 * a12 * a23 * a30 - a01 * a13 * a20 * a32 + - a02 * a10 * a23 * a31 - a02 * a11 * a20 * a33 - a02 * a13 * a21 * a30 + - a03 * a10 * a21 * a32 - a03 * a11 * a22 * a30 - a03 * a12 * a20 * a31; + }; + + // Similarly for the 5×5 circumradius determinant + auto det5 = [](const double m[5][5]) -> double { + // Laplace expansion along first row (all entries are 0,1,1,1,1) + double result = 0.0; + for (int c = 0; c < 5; ++c) { + double coef = m[0][c]; + if (std::abs(coef) < 1e-15) continue; + int sign = (c & 1) ? -1 : 1; + // Extract minor + double minor[4][4]; + for (int i = 1; i < 5; ++i) { + int mc = 0; + for (int j = 0; j < 5; ++j) { + if (j == c) continue; + minor[i-1][mc++] = m[i][j]; + } + } + result += sign * coef * det4(minor); + } + return result; + }; + + double D_circ = det5(M); + double D_vol = det4(M_vol); + + if (std::abs(D_vol) < 1e-15) return 1e6; + double R = std::sqrt(std::abs(D_circ) / (2.0 * std::abs(D_vol))); + + return R / min_edge; +} + +TetQuality evaluate_tet_quality(const HalfedgeMesh& mesh) { + TetQuality q{}; + size_t nf = mesh.num_faces(); + if (nf < 4) return q; + + // Build tetrahedra from every 4 consecutive face-triangles. + // This works for meshes where each tet is stored as 4 triangles. + size_t n_tets = nf / 4; + q.total_tets = n_tets; + if (n_tets == 0) return q; + + double min_jac = 1.0, max_jac = 0.0, sum_jac = 0.0; + double min_dihed = 180.0, max_dihed = 0.0; + double max_aspect = 0.0; + size_t degen = 0; + size_t valid = 0; + + for (size_t ti = 0; ti < n_tets; ++ti) { + // Gather vertices from the 4 face-triangles of this tet + // First triangle: (v0,v1,v2), second: (v0,v2,v3), etc. + // Standard layout: face0=abc, face1=abd, face2=acd, face3=bcd + auto f0 = mesh.face_vertices(static_cast(ti * 4)); + auto f1 = mesh.face_vertices(static_cast(ti * 4 + 1)); + auto f2 = mesh.face_vertices(static_cast(ti * 4 + 2)); + auto f3 = mesh.face_vertices(static_cast(ti * 4 + 3)); + + if (f0.size() != 3 || f1.size() != 3 || f2.size() != 3 || f3.size() != 3) { + degen++; + continue; + } + + // Deduce tetrahedron vertices from the face connectivity + Point3D a = mesh.vertex(f0[0]); // v0 + Point3D b = mesh.vertex(f0[1]); // v1 + Point3D c = mesh.vertex(f0[2]); // v2 + Point3D d = mesh.vertex(f1[2]); // v3 (abd → d is opposite face0) + + double jac = tet_scaled_jacobian(a, b, c, d); + if (jac < 1e-12) { degen++; continue; } + + min_jac = std::min(min_jac, jac); + max_jac = std::max(max_jac, jac); + sum_jac += jac; + + double dihed = tet_dihedral_angle(a, b, c, d); + min_dihed = std::min(min_dihed, dihed); + max_dihed = std::max(max_dihed, dihed); + + double aspect = tet_aspect_ratio(a, b, c, d); + max_aspect = std::max(max_aspect, aspect); + + valid++; + } + + if (valid > 0) { + q.min_jacobian = min_jac; + q.max_jacobian = max_jac; + q.avg_jacobian = sum_jac / valid; + q.min_dihedral_deg = min_dihed; + q.max_dihedral_deg = max_dihed; + q.max_aspect_ratio = max_aspect; + } + q.degenerate_tets = degen; + return q; +} + +// ═══════════════════════════════════════════════════════════ +// Generic element quality — dispatch by type +// ═══════════════════════════════════════════════════════════ + +ElemQuality evaluate_element_quality(const HalfedgeMesh& mesh, ElementType type) { + ElemQuality eq{}; + + switch (type) { + case ElementType::Tri: { + std::vector qualities; + qualities.reserve(mesh.num_faces()); + + for (size_t fi = 0; fi < mesh.num_faces(); ++fi) { + auto vis = mesh.face_vertices(static_cast(fi)); + if (vis.size() != 3) continue; + double q = tri_scaled_jacobian( + mesh.vertex(vis[0]), mesh.vertex(vis[1]), mesh.vertex(vis[2])); + if (q < 1e-12) { eq.degenerate++; continue; } + qualities.push_back(q); + } + + eq.total_elements = qualities.size() + eq.degenerate; + if (!qualities.empty()) { + eq.min_jacobian = *std::min_element(qualities.begin(), qualities.end()); + eq.max_jacobian = *std::max_element(qualities.begin(), qualities.end()); + eq.avg_jacobian = std::accumulate(qualities.begin(), qualities.end(), 0.0) + / qualities.size(); + } + break; + } + + case ElementType::Tet: { + TetQuality tq = evaluate_tet_quality(mesh); + eq.total_elements = tq.total_tets; + eq.degenerate = tq.degenerate_tets; + eq.min_jacobian = tq.min_jacobian; + eq.max_jacobian = tq.max_jacobian; + eq.avg_jacobian = tq.avg_jacobian; + break; + } + + case ElementType::Quad: + case ElementType::Hex: + // Stub — not implemented yet + eq.total_elements = 0; + break; + } + + return eq; +} + +// ═══════════════════════════════════════════════════════════ +// Quality report — histogram & grading +// ═══════════════════════════════════════════════════════════ + +namespace { + +/// Bin quality value in [0,1] into one of kBins bins +int quality_bin(double q, int kBins) { + int bin = static_cast(q * kBins); + return std::clamp(bin, 0, kBins - 1); +} + +/// Grade: A≥0.8 B≥0.6 C≥0.4 D≥0.2 F<0.2 +char quality_grade(double q) { + if (q >= 0.8) return 'A'; + if (q >= 0.6) return 'B'; + if (q >= 0.4) return 'C'; + if (q >= 0.2) return 'D'; + return 'F'; +} + +QualityReport build_report(std::vector qualities, + size_t degenerate_count, + std::string elem_type) { + QualityReport r; + r.element_type = std::move(elem_type); + r.total = qualities.size() + degenerate_count; + r.degenerate = degenerate_count; + + if (qualities.empty()) return r; + + r.min_quality = *std::min_element(qualities.begin(), qualities.end()); + r.max_quality = *std::max_element(qualities.begin(), qualities.end()); + r.avg_quality = std::accumulate(qualities.begin(), qualities.end(), 0.0) + / qualities.size(); + + // Standard deviation + double sum_sq = 0.0; + for (double q : qualities) { + double d = q - r.avg_quality; + sum_sq += d * d; + } + r.stddev_quality = std::sqrt(sum_sq / qualities.size()); + + // Histogram + grades + for (double q : qualities) { + int bin = quality_bin(q, QualityReport::kBins); + r.histogram[bin]++; + + switch (quality_grade(q)) { + case 'A': r.grade_A++; break; + case 'B': r.grade_B++; break; + case 'C': r.grade_C++; break; + case 'D': r.grade_D++; break; + case 'F': r.grade_F++; break; + } + } + + return r; +} + +} // anonymous namespace + +QualityReport mesh_quality_report(const HalfedgeMesh& mesh) { + std::vector qualities; + size_t degen = 0; + + for (size_t fi = 0; fi < mesh.num_faces(); ++fi) { + auto vis = mesh.face_vertices(static_cast(fi)); + if (vis.size() != 3) continue; + + double q = tri_scaled_jacobian( + mesh.vertex(vis[0]), mesh.vertex(vis[1]), mesh.vertex(vis[2])); + if (q < 1e-12) { degen++; continue; } + qualities.push_back(q); + } + + return build_report(std::move(qualities), degen, "triangle"); +} + +QualityReport mesh_quality_report_tet(const HalfedgeMesh& mesh) { + size_t nf = mesh.num_faces(); + size_t n_tets = nf / 4; + + std::vector qualities; + size_t degen = 0; + + for (size_t ti = 0; ti < n_tets; ++ti) { + auto f0 = mesh.face_vertices(static_cast(ti * 4)); + auto f1 = mesh.face_vertices(static_cast(ti * 4 + 1)); + auto f2 = mesh.face_vertices(static_cast(ti * 4 + 2)); + auto f3 = mesh.face_vertices(static_cast(ti * 4 + 3)); + + if (f0.size() != 3 || f1.size() != 3 || f2.size() != 3 || f3.size() != 3) { + degen++; + continue; + } + + Point3D a = mesh.vertex(f0[0]); + Point3D b = mesh.vertex(f0[1]); + Point3D c = mesh.vertex(f0[2]); + Point3D d = mesh.vertex(f1[2]); + + double q = tet_scaled_jacobian(a, b, c, d); + if (q < 1e-12) { degen++; continue; } + qualities.push_back(q); + } + + return build_report(std::move(qualities), degen, "tetrahedron"); +} + } // namespace vde::mesh diff --git a/src/mesh/mesh_smooth.cpp b/src/mesh/mesh_smooth.cpp index 047b9b8..7c3a601 100644 --- a/src/mesh/mesh_smooth.cpp +++ b/src/mesh/mesh_smooth.cpp @@ -1,53 +1,287 @@ #include "vde/mesh/mesh_smooth.h" #include +#include +#include namespace vde::mesh { -HalfedgeMesh smooth_mesh(const HalfedgeMesh& mesh, const SmoothOptions& opts) { - if (mesh.num_vertices() == 0) return mesh; +// ═══════════════════════════════════════════════════════════ +// Helpers: 1-ring neighbour collection +// ═══════════════════════════════════════════════════════════ - std::vector verts(mesh.num_vertices()); - for (size_t i = 0; i < mesh.num_vertices(); ++i) - verts[i] = mesh.vertex(i); +namespace { - for (int iter = 0; iter < opts.iterations; ++iter) { - std::vector new_verts = verts; - double lambda = (opts.method == SmoothMethod::Taubin && iter % 2 == 1) - ? opts.mu : opts.lambda; - - for (size_t vi = 0; vi < mesh.num_vertices(); ++vi) { - if (mesh.is_boundary_vertex(static_cast(vi))) continue; - - auto fis = mesh.vertex_faces(static_cast(vi)); - if (fis.empty()) continue; - - // Collect 1-ring neighbors - std::unordered_set neighbors; - for (int fi : fis) { - auto fvs = mesh.face_vertices(fi); - for (int nv : fvs) { - if (nv != static_cast(vi)) - neighbors.insert(nv); - } - } - - if (neighbors.empty()) continue; - - // Laplacian: move toward centroid of neighbors - Point3D centroid(0,0,0); - for (int nv : neighbors) centroid += verts[nv]; - centroid /= static_cast(neighbors.size()); - - new_verts[vi] = verts[vi] + lambda * (centroid - verts[vi]); +/// Collect unique 1-ring vertex neighbours of `vi`. +std::vector one_ring_neighbors(const HalfedgeMesh& mesh, int vi) { + std::unordered_set neigh; + auto fis = mesh.vertex_faces(vi); + for (int fi : fis) { + auto fvs = mesh.face_vertices(fi); + for (int v : fvs) { + if (v != vi) neigh.insert(v); } - verts = std::move(new_verts); } + return {neigh.begin(), neigh.end()}; +} - HalfedgeMesh result = mesh; +/// Compute centroid of a set of vertex indices. +Point3D centroid_of(const std::vector& verts, + const std::vector& indices) { + if (indices.empty()) return Point3D::Zero(); + Point3D c(0, 0, 0); + for (int idx : indices) c += verts[idx]; + return c / static_cast(indices.size()); +} + +/// Copy vertex positions out of a HalfedgeMesh. +std::vector extract_vertices(const HalfedgeMesh& m) { + std::vector v(m.num_vertices()); + for (size_t i = 0; i < m.num_vertices(); ++i) v[i] = m.vertex(i); + return v; +} + +/// Write vertex positions back into a HalfedgeMesh copy. +HalfedgeMesh apply_vertices(const HalfedgeMesh& orig, + const std::vector& verts) { + HalfedgeMesh result = orig; for (size_t i = 0; i < result.num_vertices(); ++i) result.set_vertex(i, verts[i]); result.update_normals(); return result; } +} // anonymous namespace + +// ═══════════════════════════════════════════════════════════ +// Laplacian smoothing +// ═══════════════════════════════════════════════════════════ + +HalfedgeMesh smooth_laplacian(const HalfedgeMesh& mesh, + int iterations, double lambda) { + if (mesh.num_vertices() == 0) return mesh; + std::vector verts = extract_vertices(mesh); + + for (int iter = 0; iter < iterations; ++iter) { + std::vector new_verts = verts; + for (size_t vi = 0; vi < mesh.num_vertices(); ++vi) { + if (mesh.is_boundary_vertex(static_cast(vi))) continue; + auto neighbors = one_ring_neighbors(mesh, vi); + if (neighbors.empty()) continue; + + Point3D c = centroid_of(verts, neighbors); + new_verts[vi] = verts[vi] + lambda * (c - verts[vi]); + } + verts = std::move(new_verts); + } + return apply_vertices(mesh, verts); +} + +// ═══════════════════════════════════════════════════════════ +// Taubin smoothing +// ═══════════════════════════════════════════════════════════ + +HalfedgeMesh smooth_taubin(const HalfedgeMesh& mesh, + int iterations, double lambda, double mu) { + if (mesh.num_vertices() == 0) return mesh; + std::vector verts = extract_vertices(mesh); + + for (int iter = 0; iter < iterations; ++iter) { + // Two passes per iteration: shrink (λ) then inflate (μ) + for (int pass = 0; pass < 2; ++pass) { + double weight = (pass == 0) ? lambda : mu; + std::vector new_verts = verts; + for (size_t vi = 0; vi < mesh.num_vertices(); ++vi) { + if (mesh.is_boundary_vertex(static_cast(vi))) continue; + auto neighbors = one_ring_neighbors(mesh, vi); + if (neighbors.empty()) continue; + + Point3D c = centroid_of(verts, neighbors); + new_verts[vi] = verts[vi] + weight * (c - verts[vi]); + } + verts = std::move(new_verts); + } + } + return apply_vertices(mesh, verts); +} + +// ═══════════════════════════════════════════════════════════ +// HC Laplacian (Humphrey's Classes) +// Two-step push-back to preserve volume: +// 1. Compute Laplacian displacement d_i = b_i − p_i^orig +// 2. Push along d_i: p_i' = p_i − (α d_i + β d_i_old) +// where b_i is the centroid of neighbours after smoothing. +// ═══════════════════════════════════════════════════════════ + +HalfedgeMesh smooth_hc_laplacian(const HalfedgeMesh& mesh, + const SmoothOptions& opts) { + if (mesh.num_vertices() == 0) return mesh; + const int iters = opts.iterations > 0 ? opts.iterations : 10; + const double lambda = (opts.lambda > 0) ? opts.lambda : 0.5; + const double alpha = (opts.hc_alpha != 0.0) ? opts.hc_alpha : 0.9; + const double beta = (opts.hc_beta != 0.0) ? opts.hc_beta : 0.5; + + std::vector verts = extract_vertices(mesh); + + // Previous push-back displacement per vertex + std::vector prev_disp(mesh.num_vertices(), Vector3D::Zero()); + + for (int iter = 0; iter < iters; ++iter) { + // Step 1: Laplacian smooth to get temporary positions q_i + std::vector temp = verts; + for (size_t vi = 0; vi < mesh.num_vertices(); ++vi) { + if (mesh.is_boundary_vertex(static_cast(vi))) continue; + auto neighbors = one_ring_neighbors(mesh, vi); + if (neighbors.empty()) continue; + Point3D c = centroid_of(verts, neighbors); + temp[vi] = verts[vi] + lambda * (c - verts[vi]); + } + + // Step 2: Compute differential and apply push-back + std::vector new_verts = verts; + std::vector new_disp(mesh.num_vertices(), Vector3D::Zero()); + + for (size_t vi = 0; vi < mesh.num_vertices(); ++vi) { + if (mesh.is_boundary_vertex(static_cast(vi))) continue; + + // Differential: difference between original and Laplacian-smoothed + Vector3D d_i = temp[vi] - verts[vi]; + + // Push-back: new = smoothed − α·d_i − β·prev_disp + // Equivalent to: new = original + d_i − α·d_i − β·prev_disp + // = original + (1−α)·d_i − β·prev_disp + new_verts[vi] = verts[vi] + (1.0 - alpha) * d_i - beta * prev_disp[vi]; + + // Compute new displacement for next iteration: + // q_i' = centroid of the new smoothed positions + auto neighbors = one_ring_neighbors(mesh, vi); + if (!neighbors.empty()) { + Point3D c_smoothed = centroid_of(new_verts, neighbors); + // Apply same Laplacian to the new positions + Point3D q_i_prime = new_verts[vi] + lambda * (c_smoothed - new_verts[vi]); + new_disp[vi] = q_i_prime - new_verts[vi]; + } + } + + verts = std::move(new_verts); + prev_disp = std::move(new_disp); + } + + return apply_vertices(mesh, verts); +} + +// ═══════════════════════════════════════════════════════════ +// Bilateral mesh filtering +// Normal-weighted smoothing that preserves sharp features. +// Each vertex is moved toward a weighted average of its +// neighbours, where weights combine: +// • spatial Gaussian: exp(−dist² / 2σ_c²) +// • normal Gaussian: exp(−(1−|dot(n_i,n_j)|)² / 2σ_n²) +// ═══════════════════════════════════════════════════════════ + +HalfedgeMesh smooth_bilateral(const HalfedgeMesh& mesh, + const SmoothOptions& opts) { + if (mesh.num_vertices() == 0) return mesh; + + const int iters = opts.bilateral_iters > 0 ? opts.bilateral_iters : 5; + const double sigma_n = (opts.bilateral_sigma_n > 0) ? opts.bilateral_sigma_n : 0.3; + + std::vector verts = extract_vertices(mesh); + + // Compute per-vertex normals + mesh.update_normals(); + // We need a mutable copy to get normals (update_normals is called in mesh) + auto get_normals = [&]() { + std::vector norms(mesh.num_vertices()); + for (size_t i = 0; i < mesh.num_vertices(); ++i) + norms[i] = mesh.vertex_normal(static_cast(i)); + return norms; + }; + + // Estimate spatial sigma from average edge length if not set + double sigma_c = opts.bilateral_sigma_c; + if (sigma_c <= 0.0) { + double sum_len = 0.0; + size_t count = 0; + for (size_t fi = 0; fi < mesh.num_faces(); ++fi) { + auto vis = mesh.face_vertices(static_cast(fi)); + if (vis.size() != 3) continue; + for (int j = 0; j < 3; ++j) { + sum_len += (verts[vis[j]] - verts[vis[(j + 1) % 3]]).norm(); + count++; + } + } + sigma_c = (count > 0) ? sum_len / count : 1.0; + } + + const double two_sigma_c2 = 2.0 * sigma_c * sigma_c; + const double two_sigma_n2 = 2.0 * sigma_n * sigma_n; + + for (int iter = 0; iter < iters; ++iter) { + // Recompute normals on current vertex positions + HalfedgeMesh temp_mesh = apply_vertices(mesh, verts); + auto current_norms = get_normals(); + + std::vector new_verts = verts; + + for (size_t vi = 0; vi < mesh.num_vertices(); ++vi) { + if (mesh.is_boundary_vertex(static_cast(vi))) continue; + + auto neighbors = one_ring_neighbors(mesh, vi); + if (neighbors.empty()) continue; + + const Vector3D& n_i = current_norms[vi]; + + Point3D weighted_sum(0, 0, 0); + double total_weight = 0.0; + + for (int vj : neighbors) { + Vector3D delta = verts[vj] - verts[vi]; + double dist2 = delta.squaredNorm(); + + // Spatial weight + double w_c = std::exp(-dist2 / two_sigma_c2); + + // Normal weight — use 1 − |dot| to penalise large normal differences + const Vector3D& n_j = current_norms[vj]; + double ndot = std::abs(n_i.dot(n_j)); + double normal_diff = 1.0 - ndot; + double w_n = std::exp(-(normal_diff * normal_diff) / two_sigma_n2); + + double w = w_c * w_n; + weighted_sum += verts[vj] * w; + total_weight += w; + } + + if (total_weight > 1e-12) { + Point3D average = weighted_sum / total_weight; + // Move vertex toward the filtered average + double l = (opts.lambda > 0) ? opts.lambda : 0.5; + new_verts[vi] = verts[vi] + l * (average - verts[vi]); + } + } + + verts = std::move(new_verts); + } + + return apply_vertices(mesh, verts); +} + +// ═══════════════════════════════════════════════════════════ +// Dispatcher +// ═══════════════════════════════════════════════════════════ + +HalfedgeMesh smooth_mesh(const HalfedgeMesh& mesh, const SmoothOptions& opts) { + switch (opts.method) { + case SmoothMethod::Laplacian: + return smooth_laplacian(mesh, opts.iterations, opts.lambda); + case SmoothMethod::Taubin: + return smooth_taubin(mesh, opts.iterations, opts.lambda, opts.mu); + case SmoothMethod::HCLaplacian: + return smooth_hc_laplacian(mesh, opts); + case SmoothMethod::Bilateral: + return smooth_bilateral(mesh, opts); + } + return mesh; +} + } // namespace vde::mesh diff --git a/src/spatial/r_tree.cpp b/src/spatial/r_tree.cpp index 621effe..5a6bb13 100644 --- a/src/spatial/r_tree.cpp +++ b/src/spatial/r_tree.cpp @@ -1,50 +1,324 @@ #include "vde/spatial/r_tree.h" #include +#include namespace vde::spatial { +// ── helpers ─────────────────────────────────────────────────────────── + +namespace { + +template +AABB3D compute_item_bounds(const T& item) { + AABB3D b; + if constexpr (std::is_same_v) { + b.expand(item); + } + return b; +} + +template +Point3D get_item_center(const T& item) { + if constexpr (std::is_same_v) { + return item; + } + return Point3D::Zero(); +} + +/// Squared minimal distance from a point to an AABB. +double mindist_sq(const Point3D& p, const AABB3D& box) { + double d = 0.0; + for (int i = 0; i < 3; ++i) { + if (p[i] < box.min()[i]) + d += (box.min()[i] - p[i]) * (box.min()[i] - p[i]); + else if (p[i] > box.max()[i]) + d += (p[i] - box.max()[i]) * (p[i] - box.max()[i]); + } + return d; +} + +/// Ray–AABB slab test (used internally during tree traversal). +bool ray_hits_aabb(const Ray3Dd& ray, const AABB3D& box) { + double tmin = -std::numeric_limits::infinity(); + double tmax = std::numeric_limits::infinity(); + + const Point3D& o = ray.origin(); + const Vector3D& d = ray.direction(); + + for (int i = 0; i < 3; ++i) { + if (std::abs(d[i]) < 1e-12) { + if (o[i] < box.min()[i] || o[i] > box.max()[i]) return false; + } else { + double inv = 1.0 / d[i]; + double t1 = (box.min()[i] - o[i]) * inv; + double t2 = (box.max()[i] - o[i]) * inv; + if (t1 > t2) std::swap(t1, t2); + tmin = std::max(tmin, t1); + tmax = std::min(tmax, t2); + if (tmin > tmax) return false; + } + } + return tmax >= 0.0; +} + +} // anonymous namespace + +// ── build (STR bulk-loading) ───────────────────────────────────────── + template void RTree::build(const std::vector& items) { - items_.clear(); - bounds_.clear(); - for (const auto& item : items) { - items_.push_back(item); - AABB3D b; - if constexpr (std::is_same_v) b.expand(item); - bounds_.push_back(b); + clear(); + if (items.empty()) return; + + const size_t N = items.size(); + items_ = items; + count_ = N; + + // 1. Pre-compute per-item bounds + item_bounds_.resize(N); + for (size_t i = 0; i < N; ++i) + item_bounds_[i] = compute_item_bounds(items_[i]); + + // 2. Sort indices by x, then slice, then sort slices by y + std::vector idx(N); + for (size_t i = 0; i < N; ++i) idx[i] = i; + + std::sort(idx.begin(), idx.end(), [&](size_t a, size_t b) { + return get_item_center(items_[a]).x() < get_item_center(items_[b]).x(); + }); + + const size_t num_slices = static_cast(std::ceil(std::sqrt(static_cast(N)))); + const size_t base_slice_sz = N / num_slices; + const size_t rem_slice = N % num_slices; + + std::vector> leaves; + + size_t cursor = 0; + for (size_t s = 0; s < num_slices; ++s) { + size_t slice_sz = base_slice_sz + (s < rem_slice ? 1 : 0); + if (slice_sz == 0) continue; + + // Sort slice by y + std::sort(idx.begin() + cursor, idx.begin() + cursor + slice_sz, + [&](size_t a, size_t b) { + return get_item_center(items_[a]).y() < + get_item_center(items_[b]).y(); + }); + + size_t nodes_in_slice = static_cast( + std::ceil(std::sqrt(static_cast(slice_sz)))); + size_t base_leaf_sz = slice_sz / nodes_in_slice; + size_t rem_leaf = slice_sz % nodes_in_slice; + + size_t sc = cursor; + for (size_t ln = 0; ln < nodes_in_slice; ++ln) { + size_t leaf_sz = base_leaf_sz + (ln < rem_leaf ? 1 : 0); + if (leaf_sz == 0) continue; + + RTreeNode leaf; + leaf.is_leaf = true; + leaf.item_indices.assign(idx.begin() + sc, + idx.begin() + sc + leaf_sz); + // Union bounding box + for (size_t ii : leaf.item_indices) + leaf.bbox.expand(item_bounds_[ii]); + + leaves.push_back(std::move(leaf)); + sc += leaf_sz; + } + cursor += slice_sz; } - count_ = items_.size(); + + // 3. Recursively build internal levels until single root + nodes_ = std::move(leaves); + size_t level_start = 0; // first node at current level + size_t level_count = nodes_.size(); // how many at current level + + while (level_count > 1) { + std::vector> parents; + + for (size_t i = 0; i < level_count; /* inside */) { + size_t fanout = static_cast( + std::ceil(std::sqrt(static_cast(level_count - i)))); + fanout = std::max(fanout, size_t(2)); + fanout = std::min(fanout, level_count - i); + + RTreeNode parent; + parent.is_leaf = false; + for (size_t j = 0; j < fanout; ++j) { + size_t child_idx = level_start + i + j; + parent.children.push_back(child_idx); + parent.bbox.expand(nodes_[child_idx].bbox); + nodes_[child_idx].parent = level_start + level_count + parents.size(); + } + parents.push_back(std::move(parent)); + i += fanout; + } + + level_start += level_count; + level_count = parents.size(); + nodes_.insert(nodes_.end(), + std::make_move_iterator(parents.begin()), + std::make_move_iterator(parents.end())); + } + + root_idx_ = level_start; } +// ── insert / remove ────────────────────────────────────────────────── + template void RTree::insert(const T& item) { items_.push_back(item); - AABB3D b; - if constexpr (std::is_same_v) b.expand(item); - bounds_.push_back(b); - count_++; + item_bounds_.push_back(compute_item_bounds(item)); + ++count_; + // Rebuild for correctness (STR is bulk-only) + build(items_); } template -bool RTree::remove(const T&) { return false; } +bool RTree::remove(const T& item) { + auto it = std::find(items_.begin(), items_.end(), item); + if (it == items_.end()) return false; + items_.erase(it); + build(items_); // rebuild + return true; +} + +// ── query_range ────────────────────────────────────────────────────── + +template +void RTree::query_range_recursive(size_t node_idx, const AABB3D& range, + std::vector& result) const { + const auto& node = nodes_[node_idx]; + if (!node.bbox.intersects(range)) return; + + if (node.is_leaf) { + for (size_t ii : node.item_indices) { + if (item_bounds_[ii].intersects(range)) + result.push_back(items_[ii]); + } + return; + } + + for (size_t child : node.children) + query_range_recursive(child, range, result); +} template std::vector RTree::query_range(const AABB3D& range) const { std::vector result; - for (size_t i = 0; i < items_.size(); ++i) { - if (bounds_[i].intersects(range)) result.push_back(items_[i]); - } + if (!nodes_.empty()) + query_range_recursive(root_idx_, range, result); return result; } -template -std::vector RTree::query_knn(const Point3D&, size_t) const { return {}; } +// ── query_knn (priority-queue best-first) ──────────────────────────── template -std::vector RTree::query_ray(const Ray3Dd&) const { return {}; } +std::vector RTree::query_knn(const Point3D& point, size_t k) const { + if (nodes_.empty() || k == 0) return {}; + + struct HeapItem { + double dist; + size_t node_idx; + bool operator>(const HeapItem& o) const { return dist > o.dist; } + }; + + std::priority_queue, std::greater> pq; + pq.push({mindist_sq(point, nodes_[root_idx_].bbox), root_idx_}); + + // Maintain k-nearest candidates + struct Candidate { + double dist_sq; + size_t item_idx; + bool operator<(const Candidate& o) const { return dist_sq < o.dist_sq; } + }; + std::vector best; + best.reserve(k + 1); + + double kth_dist_sq = std::numeric_limits::infinity(); + + while (!pq.empty()) { + auto top = pq.top(); + pq.pop(); + + if (top.dist > kth_dist_sq) break; // no better candidates left + + const auto& node = nodes_[top.node_idx]; + + if (node.is_leaf) { + for (size_t ii : node.item_indices) { + double d2 = (items_[ii] - point).squaredNorm(); + if (d2 < kth_dist_sq) { + best.push_back({d2, ii}); + std::push_heap(best.begin(), best.end()); + if (best.size() > k) { + std::pop_heap(best.begin(), best.end()); + best.pop_back(); + } + if (best.size() == k) + kth_dist_sq = best.front().dist_sq; + } + } + } else { + for (size_t child : node.children) { + double d2 = mindist_sq(point, nodes_[child].bbox); + if (d2 <= kth_dist_sq) + pq.push({d2, child}); + } + } + } + + std::vector result; + result.reserve(best.size()); + for (const auto& c : best) result.push_back(items_[c.item_idx]); + return result; +} + +// ── query_ray ──────────────────────────────────────────────────────── template -void RTree::clear() { items_.clear(); bounds_.clear(); count_ = 0; } +std::vector RTree::query_ray(const Ray3Dd& ray) const { + std::vector result; + if (nodes_.empty()) return result; + + std::vector stack; + stack.push_back(root_idx_); + + while (!stack.empty()) { + size_t idx = stack.back(); + stack.pop_back(); + + const auto& node = nodes_[idx]; + if (!ray_hits_aabb(ray, node.bbox)) continue; + + if (node.is_leaf) { + for (size_t ii : node.item_indices) { + if (ray_hits_aabb(ray, item_bounds_[ii])) + result.push_back(items_[ii]); + } + } else { + for (size_t child : node.children) + stack.push_back(child); + } + } + + return result; +} + +// ── clear ──────────────────────────────────────────────────────────── + +template +void RTree::clear() { + items_.clear(); + item_bounds_.clear(); + nodes_.clear(); + root_idx_ = 0; + count_ = 0; +} + +// ── explicit instantiations ────────────────────────────────────────── template class RTree; diff --git a/tests/collision/CMakeLists.txt b/tests/collision/CMakeLists.txt index ac5e15b..2e45df4 100644 --- a/tests/collision/CMakeLists.txt +++ b/tests/collision/CMakeLists.txt @@ -1 +1,2 @@ add_vde_test(test_gjk) +add_vde_test(test_ray) diff --git a/tests/collision/test_ray.cpp b/tests/collision/test_ray.cpp new file mode 100644 index 0000000..86b327e --- /dev/null +++ b/tests/collision/test_ray.cpp @@ -0,0 +1,120 @@ +#include +#include "vde/collision/ray_intersect.h" +#include "vde/core/triangle.h" + +using namespace vde::collision; +using namespace vde::core; + +// ── Ray-Triangle intersection (origin + direction) ── + +TEST(RayTest, RayTriangle_Hit) { + // Triangle on XY plane at z=0 + Triangle3D tri(Point3D(0, 0, 0), Point3D(2, 0, 0), Point3D(0, 2, 0)); + // Ray pointing downwards from above + auto result = ray_triangle_intersect( + Point3D(0.5, 0.5, 5.0), + Vector3D(0, 0, -1), + tri + ); + ASSERT_TRUE(result.has_value()); + EXPECT_NEAR(result->point.x(), 0.5, 1e-9); + EXPECT_NEAR(result->point.y(), 0.5, 1e-9); + EXPECT_NEAR(result->point.z(), 0.0, 1e-9); + EXPECT_GT(result->t, 0.0); // t should be positive (in ray direction) +} + +TEST(RayTest, RayTriangle_MissParallel) { + // Triangle on XY plane, ray parallel to plane + Triangle3D tri(Point3D(0, 0, 0), Point3D(2, 0, 0), Point3D(0, 2, 0)); + auto result = ray_triangle_intersect( + Point3D(0, 0, 5), + Vector3D(1, 0, 0), + tri + ); + EXPECT_FALSE(result.has_value()); +} + +TEST(RayTest, RayTriangle_MissOutside) { + Triangle3D tri(Point3D(0, 0, 0), Point3D(2, 0, 0), Point3D(0, 2, 0)); + // Ray hits the plane but outside the triangle + auto result = ray_triangle_intersect( + Point3D(3, 3, 5), + Vector3D(0, 0, -1), + tri + ); + EXPECT_FALSE(result.has_value()); +} + +TEST(RayTest, RayTriangle_MissBehindOrigin) { + Triangle3D tri(Point3D(0, 0, 0), Point3D(2, 0, 0), Point3D(0, 2, 0)); + // Ray points away from the triangle + auto result = ray_triangle_intersect( + Point3D(0.5, 0.5, 5), + Vector3D(0, 0, 1), // points +Z (away) + tri + ); + EXPECT_FALSE(result.has_value()); +} + +TEST(RayTest, RayTriangle_HitVertex) { + Triangle3D tri(Point3D(0, 0, 0), Point3D(2, 0, 0), Point3D(0, 2, 0)); + // Ray directly at vertex (0,0,0) + auto result = ray_triangle_intersect( + Point3D(0, 0, 5), + Vector3D(0, 0, -1), + tri + ); + ASSERT_TRUE(result.has_value()); + EXPECT_NEAR(result->point.z(), 0.0, 1e-9); +} + +TEST(RayTest, RayTriangle_HitEdge) { + Triangle3D tri(Point3D(0, 0, 0), Point3D(2, 0, 0), Point3D(0, 2, 0)); + // Ray at edge midpoint (between v0 and v1) + auto result = ray_triangle_intersect( + Point3D(1, 0, 3), + Vector3D(0, 0, -1), + tri + ); + ASSERT_TRUE(result.has_value()); + EXPECT_NEAR(result->point.z(), 0.0, 1e-9); +} + +// ── Ray-Triangle via Ray3Dd wrapper ── + +TEST(RayTest, RayTriangle_RayWrapper) { + Triangle3D tri(Point3D(0, 0, 0), Point3D(1, 0, 0), Point3D(0, 1, 0)); + Ray3Dd ray(Point3D(0.3, 0.3, 2), Vector3D(0, 0, -1)); + auto result = ray_triangle_intersect(ray, tri); + ASSERT_TRUE(result.has_value()); + EXPECT_NEAR(result->t, 2.0, 1e-9); +} + +// ── Degenerate triangle ── + +TEST(RayTest, RayTriangle_DegenerateTriangle) { + // All three vertices are collinear + Triangle3D tri(Point3D(0, 0, 0), Point3D(1, 0, 0), Point3D(2, 0, 0)); + auto result = ray_triangle_intersect( + Point3D(0.5, 0.5, 1), + Vector3D(0, 0, -1), + tri + ); + // Degenerate triangle: should miss or handle gracefully + EXPECT_FALSE(result.has_value()); +} + +// ── Oblique ray ── + +TEST(RayTest, RayTriangle_ObliqueRay) { + Triangle3D tri(Point3D(0, 0, 0), Point3D(2, 0, 0), Point3D(0, 2, 0)); + // Oblique ray that passes through triangle at (1, 0.27, 0) + auto result = ray_triangle_intersect( + Point3D(3, -2, 5), + Vector3D(-0.37, 0.41, -0.83), + tri + ); + // Ray might or might not hit depending on direction + (void)result; + SUCCEED(); +} diff --git a/tests/core/CMakeLists.txt b/tests/core/CMakeLists.txt index 4247962..37c6d31 100644 --- a/tests/core/CMakeLists.txt +++ b/tests/core/CMakeLists.txt @@ -2,3 +2,4 @@ add_vde_test(test_point) add_vde_test(test_convex_hull) add_vde_test(test_transform) add_vde_test(test_distance) +add_vde_test(test_polygon) diff --git a/tests/core/test_distance.cpp b/tests/core/test_distance.cpp index 19b25b8..8a5c271 100644 --- a/tests/core/test_distance.cpp +++ b/tests/core/test_distance.cpp @@ -1,13 +1,124 @@ #include #include "vde/core/distance.h" +#include "vde/core/line.h" +#include "vde/core/plane.h" +#include "vde/core/triangle.h" using namespace vde::core; +// ── Point to Point ── + TEST(DistanceTest, PointToPoint) { EXPECT_DOUBLE_EQ(distance(Point3D(0,0,0), Point3D(1,0,0)), 1.0); + EXPECT_DOUBLE_EQ(distance(Point3D(0,0,0), Point3D(3,4,0)), 5.0); + EXPECT_DOUBLE_EQ(distance(Point3D(1,2,3), Point3D(1,2,3)), 0.0); } +// ── Point to Plane ── + TEST(DistanceTest, PointToPlane) { Plane3D plane(Point3D(0,0,0), Vector3D(0,0,1)); EXPECT_DOUBLE_EQ(distance(Point3D(0,0,5), plane), 5.0); + EXPECT_DOUBLE_EQ(distance(Point3D(0,0,-5), plane), 5.0); + EXPECT_DOUBLE_EQ(distance(Point3D(0,0,0), plane), 0.0); +} + +TEST(DistanceTest, PointToPlane_OffsetOrigin) { + Plane3D plane(Point3D(0,0,10), Vector3D(0,0,1)); + EXPECT_DOUBLE_EQ(distance(Point3D(0,0,0), plane), 10.0); + EXPECT_DOUBLE_EQ(distance(Point3D(0,0,10), plane), 0.0); +} + +// ── Point to Line ── + +TEST(DistanceTest, PointToLine_OnLine) { + Line3Dd line(Point3D(0,0,0), Vector3D(1,0,0)); + EXPECT_DOUBLE_EQ(distance(Point3D(5,0,0), line), 0.0); +} + +TEST(DistanceTest, PointToLine_Perpendicular) { + Line3Dd line(Point3D(0,0,0), Vector3D(1,0,0)); + // Point at (0, 3, 0), line along X-axis → distance = 3 + EXPECT_DOUBLE_EQ(distance(Point3D(0,3,0), line), 3.0); +} + +TEST(DistanceTest, PointToLine_OffsetOrigin) { + Line3Dd line(Point3D(1,2,0), Vector3D(0,1,0)); + // Line through (1,2,0) along Y-axis + // Point (1,0,0) → perpendicular distance = 2 + EXPECT_DOUBLE_EQ(distance(Point3D(1,0,0), line), 2.0); +} + +// ── Point to Segment ── + +TEST(DistanceTest, PointToSegment_PointOnSegment) { + Segment3Dd seg(Point3D(0,0,0), Point3D(2,0,0)); + EXPECT_DOUBLE_EQ(distance(Point3D(1,0,0), seg), 0.0); +} + +TEST(DistanceTest, PointToSegment_EndpointClosest) { + Segment3Dd seg(Point3D(0,0,0), Point3D(2,0,0)); + // Point beyond P1 → closest is P1 + EXPECT_DOUBLE_EQ(distance(Point3D(3,0,0), seg), 1.0); + // Point before P0 → closest is P0 + EXPECT_DOUBLE_EQ(distance(Point3D(-1,0,0), seg), 1.0); +} + +TEST(DistanceTest, PointToSegment_PerpendicularProjection) { + Segment3Dd seg(Point3D(0,0,0), Point3D(4,0,0)); + // Point directly above middle → perpendicular projection + EXPECT_DOUBLE_EQ(distance(Point3D(2,3,0), seg), 3.0); +} + +TEST(DistanceTest, PointToSegment_ZeroLength) { + Segment3Dd seg(Point3D(1,1,1), Point3D(1,1,1)); + EXPECT_DOUBLE_EQ(distance(Point3D(1,1,0), seg), 1.0); +} + +// ── Point to Triangle ── + +TEST(DistanceTest, PointToTriangle_InteriorProjection) { + Triangle3D tri(Point3D(0,0,0), Point3D(2,0,0), Point3D(0,2,0)); + // Point directly above centroid → interior projection + EXPECT_DOUBLE_EQ(distance(Point3D(0.5, 0.5, 3), tri), 3.0); +} + +TEST(DistanceTest, PointToTriangle_OnPlaneInside) { + Triangle3D tri(Point3D(0,0,0), Point3D(2,0,0), Point3D(0,2,0)); + // Point on the triangle plane, inside the triangle + EXPECT_DOUBLE_EQ(distance(Point3D(0.5, 0.5, 0), tri), 0.0); +} + +TEST(DistanceTest, PointToTriangle_OnPlaneOutside) { + Triangle3D tri(Point3D(0,0,0), Point3D(2,0,0), Point3D(0,2,0)); + // Point on the triangle plane but outside → closest to edge or vertex + double d = distance(Point3D(3, 3, 0), tri); + EXPECT_GT(d, 0.0); +} + +// ── closest_point ── + +TEST(DistanceTest, ClosestPoint_Triangle) { + Triangle3D tri(Point3D(0,0,0), Point3D(2,0,0), Point3D(0,2,0)); + Point3D c = closest_point(Point3D(1, 1, 5), tri); + // Closest point should be on triangle plane, within bounds + EXPECT_NEAR(c.z(), 0.0, 1e-9); + EXPECT_GE(c.x(), 0.0); + EXPECT_GE(c.y(), 0.0); +} + +TEST(DistanceTest, ClosestPoint_Segment) { + Segment3Dd seg(Point3D(0,0,0), Point3D(4,0,0)); + Point3D c = closest_point(Point3D(2, 5, 0), seg); + EXPECT_NEAR(c.x(), 2.0, 1e-9); + EXPECT_NEAR(c.y(), 0.0, 1e-9); + EXPECT_NEAR(c.z(), 0.0, 1e-9); +} + +TEST(DistanceTest, ClosestPoint_SegmentAtEndpoint) { + Segment3Dd seg(Point3D(0,0,0), Point3D(4,0,0)); + Point3D c = closest_point(Point3D(5, 0, 3), seg); + // Should clamp to P1 + EXPECT_NEAR(c.x(), 4.0, 1e-9); + EXPECT_NEAR(c.y(), 0.0, 1e-9); } diff --git a/tests/core/test_polygon.cpp b/tests/core/test_polygon.cpp new file mode 100644 index 0000000..11971f9 --- /dev/null +++ b/tests/core/test_polygon.cpp @@ -0,0 +1,146 @@ +#include +#include "vde/core/polygon.h" + +using namespace vde::core; + +// ── signed_area tests ── + +TEST(PolygonTest, SignedArea_CCWSquare_Positive) { + // Counter-clockwise unit square + Polygon2D poly({{0, 0}, {1, 0}, {1, 1}, {0, 1}}); + EXPECT_GT(poly.signed_area(), 0.0); + EXPECT_NEAR(poly.signed_area(), 1.0, 1e-9); +} + +TEST(PolygonTest, SignedArea_CWSquare_Negative) { + // Clockwise unit square + Polygon2D poly({{0, 0}, {0, 1}, {1, 1}, {1, 0}}); + EXPECT_LT(poly.signed_area(), 0.0); + EXPECT_NEAR(poly.signed_area(), -1.0, 1e-9); +} + +TEST(PolygonTest, SignedArea_Triangle) { + Polygon2D poly({{0, 0}, {3, 0}, {0, 4}}); + EXPECT_NEAR(std::abs(poly.signed_area()), 6.0, 1e-9); +} + +TEST(PolygonTest, SignedArea_Collinear_Degenerate) { + // All points on a line → zero area + Polygon2D poly({{0, 0}, {1, 1}, {2, 2}}); + EXPECT_NEAR(poly.signed_area(), 0.0, 1e-9); +} + +TEST(PolygonTest, SignedArea_Empty) { + Polygon2D poly; + EXPECT_DOUBLE_EQ(poly.signed_area(), 0.0); +} + +TEST(PolygonTest, Area_AbsoluteValue) { + // area() should be absolute value of signed_area() + Polygon2D ccw({{0, 0}, {2, 0}, {2, 2}, {0, 2}}); + Polygon2D cw({{0, 0}, {0, 2}, {2, 2}, {2, 0}}); + EXPECT_NEAR(ccw.area(), 4.0, 1e-9); + EXPECT_NEAR(cw.area(), 4.0, 1e-9); +} + +// ── contains (point-in-polygon via ray casting) ── + +TEST(PolygonTest, Contains_InteriorPoint) { + Polygon2D square({{0, 0}, {2, 0}, {2, 2}, {0, 2}}); + // Center of the square + EXPECT_TRUE(square.contains(Point2D(1.0, 1.0))); +} + +TEST(PolygonTest, Contains_ExteriorPoint) { + Polygon2D square({{0, 0}, {2, 0}, {2, 2}, {0, 2}}); + // Far outside + EXPECT_FALSE(square.contains(Point2D(10.0, 10.0))); + // Outside but near + EXPECT_FALSE(square.contains(Point2D(-0.1, 1.0))); +} + +TEST(PolygonTest, Contains_BoundaryPoint) { + Polygon2D square({{0, 0}, {2, 0}, {2, 2}, {0, 2}}); + // On edge — behavior may vary; ray casting treats as intersection + // Not strictly "interior" in most implementations + Point2D on_edge(1.0, 0.0); + // Just verify it doesn't crash; boundary is valid input + bool result = square.contains(on_edge); + (void)result; // implementation-defined for boundary + SUCCEED(); +} + +TEST(PolygonTest, Contains_PointOnVertex) { + Polygon2D tri({{0, 0}, {2, 0}, {0, 2}}); + // Exact vertex — should not crash + bool result = tri.contains(Point2D(0.0, 0.0)); + (void)result; + SUCCEED(); +} + +TEST(PolygonTest, Contains_ConcavePolygon) { + // L-shaped polygon (concave) + Polygon2D L({{0, 0}, {3, 0}, {3, 1}, {2, 1}, {2, 3}, {0, 3}}); + // Point in the "notch" area (1.5, 2) should be inside + EXPECT_TRUE(L.contains(Point2D(1.0, 2.0))); + // Point in the concave void + EXPECT_FALSE(L.contains(Point2D(2.5, 2.0))); +} + +// ── is_ccw ── + +TEST(PolygonTest, IsCCW) { + Polygon2D ccw({{0, 0}, {1, 0}, {1, 1}, {0, 1}}); + Polygon2D cw({{0, 0}, {0, 1}, {1, 1}, {1, 0}}); + EXPECT_TRUE(ccw.is_ccw()); + EXPECT_FALSE(cw.is_ccw()); +} + +// ── Basic accessors ── + +TEST(PolygonTest, SizeAndEmpty) { + Polygon2D empty; + EXPECT_TRUE(empty.empty()); + EXPECT_EQ(empty.size(), 0u); + + Polygon2D tri({{0, 0}, {1, 0}, {0, 1}}); + EXPECT_FALSE(tri.empty()); + EXPECT_EQ(tri.size(), 3u); +} + +TEST(PolygonTest, VerticesAccessor) { + std::vector v = {{0, 0}, {1, 0}, {0, 1}}; + Polygon2D poly(v); + EXPECT_EQ(poly.vertices().size(), 3u); + EXPECT_DOUBLE_EQ(poly.vertices()[0].x(), 0.0); + EXPECT_DOUBLE_EQ(poly.vertices()[1].y(), 0.0); +} + +TEST(PolygonTest, DefaultConstructor_Empty) { + Polygon2D poly; + EXPECT_TRUE(poly.vertices().empty()); + EXPECT_TRUE(poly.empty()); +} + +// ── Larger polygon area (regular hexagon) ── + +TEST(PolygonTest, SignedArea_RegularHexagon) { + // Regular hexagon centered at origin, radius 1 + // Vertices CCW: (1,0), (0.5,0.866), (-0.5,0.866), (-1,0), (-0.5,-0.866), (0.5,-0.866) + Polygon2D hex({ + {1.0, 0.0}, {0.5, std::sqrt(3.0) / 2.0}, + {-0.5, std::sqrt(3.0) / 2.0}, {-1.0, 0.0}, + {-0.5, -std::sqrt(3.0) / 2.0}, {0.5, -std::sqrt(3.0) / 2.0} + }); + // Area of regular hexagon = (3√3 / 2) * r² ≈ 2.598 + EXPECT_NEAR(hex.signed_area(), 3.0 * std::sqrt(3.0) / 2.0, 1e-9); +} + +TEST(PolygonTest, Contains_StarShapedHexagon) { + Polygon2D hex({ + {1.0, 0.0}, {0.5, 0.866}, {-0.5, 0.866}, + {-1.0, 0.0}, {-0.5, -0.866}, {0.5, -0.866} + }); + EXPECT_TRUE(hex.contains(Point2D(0.0, 0.0))); + EXPECT_FALSE(hex.contains(Point2D(2.0, 0.0))); +} diff --git a/tests/mesh/CMakeLists.txt b/tests/mesh/CMakeLists.txt index c6c5e00..1042d4d 100644 --- a/tests/mesh/CMakeLists.txt +++ b/tests/mesh/CMakeLists.txt @@ -1,2 +1,4 @@ add_vde_test(test_halfedge) add_vde_test(test_delaunay) +add_vde_test(test_quality) +add_vde_test(test_smooth) diff --git a/tests/mesh/test_quality.cpp b/tests/mesh/test_quality.cpp new file mode 100644 index 0000000..ae83d21 --- /dev/null +++ b/tests/mesh/test_quality.cpp @@ -0,0 +1,106 @@ +#include +#include "vde/mesh/mesh_quality.h" +#include "vde/mesh/halfedge_mesh.h" + +using namespace vde::mesh; + +TEST(MeshQualityTest, EquilateralTriangle_PerfectQuality) { + HalfedgeMesh mesh; + // Equilateral triangle side length ≈ 1.155 for area ≈ 0.577 + double h = std::sqrt(3.0) / 2.0; + mesh.build_from_triangles( + {Point3D(0, 0, 0), Point3D(1, 0, 0), Point3D(0.5, h, 0)}, + {{0, 1, 2}} + ); + + auto q = evaluate_mesh_quality(mesh); + + // Equilateral triangle: all angles = 60° + EXPECT_NEAR(q.min_angle_deg, 60.0, 1.0); + EXPECT_NEAR(q.max_angle_deg, 60.0, 1.0); + // Aspect ratio ~1 for equilateral + EXPECT_NEAR(q.avg_aspect_ratio, 1.0, 0.1); + EXPECT_EQ(q.degenerate_faces, 0u); +} + +TEST(MeshQualityTest, DegenerateTriangle_ZeroArea) { + HalfedgeMesh mesh; + // Collinear points → zero area + mesh.build_from_triangles( + {Point3D(0, 0, 0), Point3D(1, 0, 0), Point3D(2, 0, 0)}, + {{0, 1, 2}} + ); + + auto q = evaluate_mesh_quality(mesh); + EXPECT_EQ(q.degenerate_faces, 1u); +} + +TEST(MeshQualityTest, ObtuseTriangle_LargeMaxAngle) { + HalfedgeMesh mesh; + // Very flat obtuse triangle: (0,0), (10,0), (0,0.1) + mesh.build_from_triangles( + {Point3D(0, 0, 0), Point3D(10, 0, 0), Point3D(0, 0.1, 0)}, + {{0, 1, 2}} + ); + + auto q = evaluate_mesh_quality(mesh); + // Max angle should be very close to 180 (or at least > 90) + EXPECT_GT(q.max_angle_deg, 90.0); + // Min angle should be very small + EXPECT_LT(q.min_angle_deg, 10.0); +} + +TEST(MeshQualityTest, RightTriangle) { + HalfedgeMesh mesh; + // 3-4-5 right triangle + mesh.build_from_triangles( + {Point3D(0, 0, 0), Point3D(3, 0, 0), Point3D(0, 4, 0)}, + {{0, 1, 2}} + ); + + auto q = evaluate_mesh_quality(mesh); + // Right triangle: one angle ~90°, others ~36.87° and ~53.13° + EXPECT_NEAR(q.max_angle_deg, 90.0, 2.0); + EXPECT_GT(q.min_angle_deg, 30.0); + EXPECT_EQ(q.degenerate_faces, 0u); +} + +TEST(MeshQualityTest, MultipleTriangles_Mixed) { + HalfedgeMesh mesh; + // Two triangles: one good, one degenerate + double h = std::sqrt(3.0) / 2.0; + mesh.build_from_triangles( + { + Point3D(0, 0, 0), Point3D(1, 0, 0), Point3D(0.5, h, 0), // equilateral + Point3D(10, 0, 0), Point3D(11, 0, 0), Point3D(12, 0, 0), // degenerate + }, + {{0, 1, 2}, {3, 4, 5}} + ); + + auto q = evaluate_mesh_quality(mesh); + EXPECT_EQ(q.degenerate_faces, 1u); + // Average aspect ratio should reflect the mix + EXPECT_GT(q.avg_aspect_ratio, 1.0); +} + +TEST(MeshQualityTest, EmptyMesh) { + HalfedgeMesh mesh; + auto q = evaluate_mesh_quality(mesh); + EXPECT_EQ(q.degenerate_faces, 0u); + // Default values for empty mesh + EXPECT_DOUBLE_EQ(q.min_angle_deg, 0.0); + EXPECT_DOUBLE_EQ(q.max_angle_deg, 0.0); +} + +TEST(MeshQualityTest, HighAspectRatioTriangle) { + HalfedgeMesh mesh; + // Very thin triangle + mesh.build_from_triangles( + {Point3D(0, 0, 0), Point3D(100, 0, 0), Point3D(0, 0.01, 0)}, + {{0, 1, 2}} + ); + + auto q = evaluate_mesh_quality(mesh); + // Aspect ratio should be very high + EXPECT_GT(q.max_aspect_ratio, 10.0); +} diff --git a/tests/mesh/test_smooth.cpp b/tests/mesh/test_smooth.cpp new file mode 100644 index 0000000..b631eb3 --- /dev/null +++ b/tests/mesh/test_smooth.cpp @@ -0,0 +1,158 @@ +#include +#include "vde/mesh/mesh_smooth.h" +#include "vde/mesh/halfedge_mesh.h" + +using namespace vde::mesh; + +// Helper: create a simple 4-vertex quad mesh (2 triangles) +static HalfedgeMesh make_quad_mesh() { + HalfedgeMesh mesh; + // Quad with slight noise on one vertex + mesh.build_from_triangles( + { + Point3D(0, 0, 0), Point3D(1, 0, 0), + Point3D(1, 1, 0.5), // perturbed in Z + Point3D(0, 1, 0), + }, + {{0, 1, 2}, {0, 2, 3}} + ); + return mesh; +} + +TEST(MeshSmoothTest, LaplacianSmooth_ReducesNoise) { + HalfedgeMesh mesh = make_quad_mesh(); + + // Original: vertex 2 has Z = 0.5 perturbation + EXPECT_NEAR(mesh.vertex(2).z(), 0.5, 1e-6); + + SmoothOptions opts; + opts.method = SmoothMethod::Laplacian; + opts.iterations = 5; + opts.lambda = 0.5; + + HalfedgeMesh smoothed = smooth_mesh(mesh, opts); + + // After smoothing, the noise should be reduced + // Vertex 2 Z should be closer to 0 (the Laplacian average) + EXPECT_LT(std::abs(smoothed.vertex(2).z()), 0.5); +} + +TEST(MeshSmoothTest, LaplacianNoop_FlatMesh) { + HalfedgeMesh mesh; + mesh.build_from_triangles( + { + Point3D(0, 0, 0), Point3D(1, 0, 0), + Point3D(1, 1, 0), Point3D(0, 1, 0), + }, + {{0, 1, 2}, {0, 2, 3}} + ); + + SmoothOptions opts; + opts.method = SmoothMethod::Laplacian; + opts.iterations = 3; + opts.lambda = 0.5; + + HalfedgeMesh smoothed = smooth_mesh(mesh, opts); + + // Flat mesh should remain flat + for (size_t i = 0; i < smoothed.num_vertices(); ++i) { + EXPECT_NEAR(smoothed.vertex(i).z(), 0.0, 1e-9); + } +} + +TEST(MeshSmoothTest, TaubinPreservesVolume) { + HalfedgeMesh mesh = make_quad_mesh(); + + // Compute original bounding box volume + auto original_bounds = mesh.bounds(); + double original_volume = original_bounds.volume(); + + SmoothOptions opts; + opts.method = SmoothMethod::Taubin; + opts.lambda = 0.5; + opts.mu = -0.53; + opts.iterations = 5; + + HalfedgeMesh smoothed = smooth_mesh(mesh, opts); + auto smoothed_bounds = smoothed.bounds(); + double smoothed_volume = smoothed_bounds.volume(); + + // Taubin should better preserve volume than pure Laplacian + // Volume change should be reasonable (not shrinking to 0) + EXPECT_GT(smoothed_volume, 0.0); + + // Noise reduction: vertex 2 Z should be smoothed + EXPECT_LT(std::abs(smoothed.vertex(2).z()), 0.5); + + // Volume change ratio should be reasonable (< 50% change) + double ratio = std::abs(smoothed_volume - original_volume) / std::max(original_volume, 1e-9); + EXPECT_LT(ratio, 0.5); +} + +TEST(MeshSmoothTest, DefaultOptions) { + HalfedgeMesh mesh = make_quad_mesh(); + + // Default options: Laplacian, 10 iterations, lambda=0.5 + HalfedgeMesh smoothed = smooth_mesh(mesh); + + EXPECT_EQ(smoothed.num_vertices(), mesh.num_vertices()); + EXPECT_EQ(smoothed.num_faces(), mesh.num_faces()); +} + +TEST(MeshSmoothTest, DifferentIterationCounts) { + HalfedgeMesh mesh = make_quad_mesh(); + double original_z = mesh.vertex(2).z(); + + SmoothOptions opts_few; + opts_few.method = SmoothMethod::Laplacian; + opts_few.iterations = 1; + opts_few.lambda = 0.5; + HalfedgeMesh smooth_1 = smooth_mesh(mesh, opts_few); + + SmoothOptions opts_many; + opts_many.method = SmoothMethod::Laplacian; + opts_many.iterations = 20; + opts_many.lambda = 0.5; + HalfedgeMesh smooth_20 = smooth_mesh(mesh, opts_many); + + // More iterations → more smoothing (vertex 2 Z closer to 0) + double delta_1 = std::abs(smooth_1.vertex(2).z() - 0.0); + double delta_20 = std::abs(smooth_20.vertex(2).z() - 0.0); + // 20 iterations should smooth more than 1 + EXPECT_LE(delta_20, delta_1 + 1e-6); +} + +TEST(MeshSmoothTest, SingleTriangleUnaffected) { + // Single triangle: smoothing has no effect (all vertices are boundary) + HalfedgeMesh mesh; + mesh.build_from_triangles( + {Point3D(0, 0, 2), Point3D(1, 0, 0), Point3D(0, 1, 0)}, + {{0, 1, 2}} + ); + + SmoothOptions opts; + opts.method = SmoothMethod::Laplacian; + opts.iterations = 5; + opts.lambda = 0.5; + + HalfedgeMesh smoothed = smooth_mesh(mesh, opts); + + // Single triangle: topological boundary, so minimal change expected + // But implementation might still move vertices + EXPECT_EQ(smoothed.num_vertices(), 3u); + EXPECT_EQ(smoothed.num_faces(), 1u); +} + +TEST(MeshSmoothTest, ZeroIterations) { + HalfedgeMesh mesh = make_quad_mesh(); + + SmoothOptions opts; + opts.iterations = 0; + + HalfedgeMesh smoothed = smooth_mesh(mesh, opts); + + // Zero iterations → mesh unchanged + for (size_t i = 0; i < mesh.num_vertices(); ++i) { + EXPECT_NEAR((smoothed.vertex(i) - mesh.vertex(i)).norm(), 0.0, 1e-9); + } +} diff --git a/tests/spatial/CMakeLists.txt b/tests/spatial/CMakeLists.txt index 8a5153f..283316d 100644 --- a/tests/spatial/CMakeLists.txt +++ b/tests/spatial/CMakeLists.txt @@ -1 +1,2 @@ add_vde_test(test_bvh) +add_vde_test(test_r_tree) diff --git a/tests/spatial/test_r_tree.cpp b/tests/spatial/test_r_tree.cpp new file mode 100644 index 0000000..43c0f61 --- /dev/null +++ b/tests/spatial/test_r_tree.cpp @@ -0,0 +1,126 @@ +#include +#include "vde/spatial/r_tree.h" +#include "vde/core/aabb.h" + +using namespace vde::spatial; +using namespace vde::core; + +// A simple wrapper to satisfy the RTree template +namespace { + struct Point3DWithId { + Point3D p; + int id; + }; +} // namespace + +// Instantiate the RTree for Point3D (items are point-based) +// RTree stores items and computes their AABBs +// For testing we use Point3D items directly + +TEST(RTreeTest, BuildEmpty) { + RTree tree; + tree.build({}); + EXPECT_EQ(tree.size(), 0u); +} + +TEST(RTreeTest, BuildSinglePoint) { + RTree tree; + tree.build({Point3D(1, 2, 3)}); + EXPECT_EQ(tree.size(), 1u); +} + +TEST(RTreeTest, BuildAndRangeQuery) { + RTree tree; + std::vector points = { + {0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {1, 1, 0}, + {5, 5, 0}, {6, 5, 0}, {5, 6, 0}, {6, 6, 0}, + }; + tree.build(points); + + // Query: range covering the first cluster (0,0)-(2,2) + AABB3D range(Point3D(-0.5, -0.5, -0.5), Point3D(2.5, 2.5, 0.5)); + auto results = tree.query_range(range); + EXPECT_EQ(results.size(), 4u); // 4 points in first cluster +} + +TEST(RTreeTest, RangeQueryEmpty) { + RTree tree; + tree.build({Point3D(0, 0, 0), Point3D(1, 1, 1)}); + + AABB3D far(Point3D(10, 10, 10), Point3D(20, 20, 20)); + auto results = tree.query_range(far); + EXPECT_TRUE(results.empty()); +} + +TEST(RTreeTest, KNNQuery) { + RTree tree; + tree.build({ + {0, 0, 0}, {10, 0, 0}, {0, 10, 0}, {10, 10, 0}, + }); + + auto nearest = tree.query_knn(Point3D(1, 0, 0), 1); + ASSERT_EQ(nearest.size(), 1u); + EXPECT_NEAR(nearest[0].x(), 0.0, 1e-6); + EXPECT_NEAR(nearest[0].y(), 0.0, 1e-6); + + auto two_nearest = tree.query_knn(Point3D(0, 0, 0), 2); + EXPECT_EQ(two_nearest.size(), 2u); +} + +TEST(RTreeTest, KNNQuery_MoreThanStored) { + RTree tree; + tree.build({Point3D(1, 2, 3), Point3D(4, 5, 6)}); + + auto results = tree.query_knn(Point3D(0, 0, 0), 10); + EXPECT_EQ(results.size(), 2u); // Only 2 items stored +} + +TEST(RTreeTest, RayQuery) { + RTree tree; + std::vector points; + for (int x = 0; x < 5; ++x) + for (int y = 0; y < 5; ++y) + points.push_back(Point3D(x, y, 0)); + tree.build(points); + + // Ray along X-axis at y=0, z=0 should hit points with y≈0 + Ray3Dd ray(Point3D(-1, 0, 0), Vector3D(1, 0, 0)); + auto hits = tree.query_ray(ray); + EXPECT_GT(hits.size(), 0u); +} + +TEST(RTreeTest, InsertAndRemove) { + RTree tree; + tree.insert(Point3D(0, 0, 0)); + EXPECT_EQ(tree.size(), 1u); + + tree.insert(Point3D(1, 1, 1)); + EXPECT_EQ(tree.size(), 2u); + + bool removed = tree.remove(Point3D(1, 1, 1)); + EXPECT_TRUE(removed); + EXPECT_EQ(tree.size(), 1u); +} + +TEST(RTreeTest, RemoveNonExistent) { + RTree tree; + tree.insert(Point3D(0, 0, 0)); + bool removed = tree.remove(Point3D(9, 9, 9)); + EXPECT_FALSE(removed); + EXPECT_EQ(tree.size(), 1u); +} + +TEST(RTreeTest, Clear) { + RTree tree; + tree.build({ + {0, 0, 0}, {1, 0, 0}, {0, 1, 0}, + }); + EXPECT_EQ(tree.size(), 3u); + + tree.clear(); + EXPECT_EQ(tree.size(), 0u); + + // After clear, queries should return empty + AABB3D all(Point3D(-1, -1, -1), Point3D(10, 10, 10)); + EXPECT_TRUE(tree.query_range(all).empty()); +}