feat: LOD mesh generation + incremental BVH

src/CMakeLists.txt

@@ -72,6 +72,7 @@ add_library(vde_mesh STATIC
     mesh/geodesic.cpp
     mesh/mesh_curvature.cpp
     mesh/marching_cubes.cpp
+    mesh/mesh_lod.cpp
 )
 target_include_directories(vde_mesh
     PUBLIC  ${CMAKE_SOURCE_DIR}/include
@@ -87,6 +88,7 @@ add_library(vde_spatial STATIC
     spatial/octree.cpp
     spatial/kd_tree.cpp
     spatial/r_tree.cpp
+    spatial/incremental_bvh.cpp
 )
 target_include_directories(vde_spatial
     PUBLIC  ${CMAKE_SOURCE_DIR}/include

src/mesh/mesh_lod.cpp

@@ -0,0 +1,55 @@
+#include "vde/mesh/mesh_lod.h"
+#include "vde/mesh/mesh_simplify.h"
+#include "vde/brep/brep.h"
+#include <algorithm>
+
+namespace vde::mesh {
+
+LODMesh generate_lod(const HalfedgeMesh& mesh, int num_levels, double reduction_ratio) {
+    LODMesh result;
+    result.name = "lod";
+    result.levels.push_back(mesh);
+    result.thresholds.push_back(0.0);  // base level: always use
+
+    HalfedgeMesh current = mesh;
+    for (int i = 1; i < num_levels; ++i) {
+        int current_faces = static_cast<int>(current.num_faces());
+        int target_faces = static_cast<int>(current_faces * reduction_ratio);
+        if (target_faces < 4) break;
+
+        SimplifyOptions opts;
+        opts.target_ratio = static_cast<double>(target_faces) / current_faces;
+        opts.max_iterations = 100;
+        opts.preserve_boundary = true;
+
+        current = simplify_mesh(current, opts);
+        result.levels.push_back(current);
+        result.thresholds.push_back(static_cast<double>(i) * 100.0);
+    }
+
+    return result;
+}
+
+LODMesh generate_brep_lod(const brep::BrepModel& body, int num_levels) {
+    auto mesh = body.to_mesh(0.1);
+    return generate_lod(mesh, num_levels);
+}
+
+int LODMesh::level_for_error(double error) const {
+    for (int i = static_cast<int>(thresholds.size()) - 1; i >= 0; --i) {
+        if (error >= thresholds[i]) return i;
+    }
+    return 0;
+}
+
+const HalfedgeMesh& LODMesh::mesh_at(int level) const {
+    if (levels.empty()) {
+        static HalfedgeMesh empty;
+        return empty;
+    }
+    if (level < 0) level = 0;
+    if (level >= static_cast<int>(levels.size())) level = static_cast<int>(levels.size()) - 1;
+    return levels[level];
+}
+
+} // namespace vde::mesh

src/spatial/incremental_bvh.cpp

@@ -0,0 +1,115 @@
+#include "vde/spatial/incremental_bvh.h"
+#include "vde/core/line.h"
+#include <cmath>
+#include <algorithm>
+
+namespace vde::spatial {
+
+using core::Point3D;
+using core::Vector3D;
+
+namespace {
+
+/// Compute AABB from triangle vertices
+AABB3D triangle_bounds(const Triangle3D& tri) {
+    double min_x = std::min({tri.v(0).x(), tri.v(1).x(), tri.v(2).x()});
+    double min_y = std::min({tri.v(0).y(), tri.v(1).y(), tri.v(2).y()});
+    double min_z = std::min({tri.v(0).z(), tri.v(1).z(), tri.v(2).z()});
+    double max_x = std::max({tri.v(0).x(), tri.v(1).x(), tri.v(2).x()});
+    double max_y = std::max({tri.v(0).y(), tri.v(1).y(), tri.v(2).y()});
+    double max_z = std::max({tri.v(0).z(), tri.v(1).z(), tri.v(2).z()});
+    return AABB3D(Point3D(min_x, min_y, min_z), Point3D(max_x, max_y, max_z));
+}
+
+/// Möller–Trumbore ray-triangle intersection
+bool mt_intersect(const Ray3Dd& ray, const Triangle3D& tri,
+                  double& t_out, Point3D& point_out) {
+    const auto& v0 = tri.v(0);
+    const auto& v1 = tri.v(1);
+    const auto& v2 = tri.v(2);
+
+    Vector3D e1 = v1 - v0;
+    Vector3D e2 = v2 - v0;
+    Vector3D pvec = ray.direction().cross(e2);
+    double det = e1.dot(pvec);
+
+    if (std::abs(det) < 1e-12) return false;
+
+    double inv_det = 1.0 / det;
+    Vector3D tvec = ray.origin() - v0;
+    double u = tvec.dot(pvec) * inv_det;
+    if (u < 0.0 || u > 1.0) return false;
+
+    Vector3D qvec = tvec.cross(e1);
+    double v = ray.direction().dot(qvec) * inv_det;
+    if (v < 0.0 || u + v > 1.0) return false;
+
+    t_out = e2.dot(qvec) * inv_det;
+    if (t_out <= 0.0) return false;
+
+    point_out = ray.origin() + ray.direction() * t_out;
+    return true;
+}
+
+} // anonymous namespace
+
+void IncrementalBVH::build(const std::vector<Triangle3D>& triangles) {
+    triangles_ = triangles;
+    bboxes_.resize(triangles.size());
+    for (size_t i = 0; i < triangles.size(); ++i) {
+        bboxes_[i] = triangle_bounds(triangles_[i]);
+    }
+    bvh_.build(triangles_);
+    dirty_ = false;
+}
+
+void IncrementalBVH::insert(const Triangle3D& tri) {
+    triangles_.push_back(tri);
+    bboxes_.push_back(triangle_bounds(tri));
+    dirty_ = true;
+}
+
+void IncrementalBVH::remove(size_t index) {
+    if (index < triangles_.size()) {
+        triangles_.erase(triangles_.begin() + static_cast<ptrdiff_t>(index));
+        bboxes_.erase(bboxes_.begin() + static_cast<ptrdiff_t>(index));
+        dirty_ = true;
+    }
+}
+
+void IncrementalBVH::update(size_t index, const Triangle3D& tri) {
+    if (index < triangles_.size()) {
+        triangles_[index] = tri;
+        bboxes_[index] = triangle_bounds(tri);
+        dirty_ = true;
+    }
+}
+
+void IncrementalBVH::rebuild_if_dirty() {
+    if (dirty_ && !triangles_.empty()) {
+        bvh_ = BVH();
+        bvh_.build(triangles_);
+        dirty_ = false;
+    }
+}
+
+std::vector<BVH::HitResult> IncrementalBVH::ray_query(const Ray3Dd& ray) const {
+    // Rebuild tree if modifications have occurred since last query
+    const_cast<IncrementalBVH*>(this)->rebuild_if_dirty();
+
+    // Use BVH to find candidate triangles
+    auto candidates = bvh_.query_ray(ray);
+
+    // Compute exact intersections for each candidate
+    std::vector<BVH::HitResult> hits;
+    for (const auto& tri : candidates) {
+        double t;
+        Point3D point;
+        if (mt_intersect(ray, tri, t, point)) {
+            hits.push_back({t, tri, point});
+        }
+    }
+    return hits;
+}
+
+} // namespace vde::spatial