diff --git a/CMakeLists.txt b/CMakeLists.txt index f3eab9b..61dcee2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -47,6 +47,7 @@ if(BUILD_TESTS) endif() # ── 示例 ────────────────────────────────────────── +# ── Python 绑定 ────────────────────────────────────option(VDE_BUILD_PYTHON "Build Python bindings (pybind11)" OFF)if(VDE_BUILD_PYTHON) add_subdirectory(python)endif() if(BUILD_EXAMPLES) add_subdirectory(examples) endif() diff --git a/include/vde/core/voronoi.h b/include/vde/core/voronoi.h new file mode 100644 index 0000000..b6caa90 --- /dev/null +++ b/include/vde/core/voronoi.h @@ -0,0 +1,18 @@ +#pragma once +#include "vde/core/point.h" +#include +#include + +namespace vde::core { + +struct VoronoiCell { + Point2D site; // Generator point + std::vector vertices; // Cell boundary vertices + std::vector> edges; // Edge indices into vertices +}; + +/// 2D Voronoi diagram from Delaunay dual graph +/// Returns cells for each input point +std::vector voronoi_2d(const std::vector& points); + +} // namespace vde::core diff --git a/include/vde/foundation/serializer.h b/include/vde/foundation/serializer.h new file mode 100644 index 0000000..ac9482a --- /dev/null +++ b/include/vde/foundation/serializer.h @@ -0,0 +1,38 @@ +#pragma once +#include "vde/mesh/halfedge_mesh.h" +#include +#include +#include + +namespace vde::foundation { + +// Versioned binary format for fast save/load +// Header: 8 bytes magic + 4 bytes version + 4 bytes flags + +constexpr uint64_t VDE_MAGIC = 0x56444547454F4D00ULL; // "VDEGEOM\0" +constexpr uint32_t VDE_FORMAT_VERSION = 1; + +struct SerializedMesh { + uint32_t vertex_count; + uint32_t face_count; + std::vector vertices; // xyz xyz ... (flat) + std::vector face_indices; // v0 v1 v2 ... + std::vector face_valence; // 3 for all triangles +}; + +class BinarySerializer { +public: + /// Serialize mesh to binary buffer + static std::vector serialize(const HalfedgeMesh& mesh); + + /// Deserialize mesh from binary buffer + static HalfedgeMesh deserialize(const std::vector& data); + + /// Write to file + static bool write_file(const std::string& path, const HalfedgeMesh& mesh); + + /// Read from file + static HalfedgeMesh read_file(const std::string& path); +}; + +} // namespace vde::foundation diff --git a/include/vde/mesh/marching_cubes.h b/include/vde/mesh/marching_cubes.h new file mode 100644 index 0000000..6f2dc93 --- /dev/null +++ b/include/vde/mesh/marching_cubes.h @@ -0,0 +1,47 @@ +#pragma once +#include "vde/core/point.h" +#include +#include +#include + +namespace vde::mesh { + +struct MCMesh { + std::vector vertices; + std::vector> triangles; +}; + +/// Marching Cubes for scalar field f(x,y,z) +/// Extracts isosurface at value = iso_level within the bounding box +MCMesh marching_cubes(const std::function& f, + double iso_level, + const Point3D& bmin, const Point3D& bmax, + int resolution); + +/// SDF sphere +inline double sdf_sphere(double x, double y, double z, double cx, double cy, double cz, double r) { + return std::sqrt((x-cx)*(x-cx) + (y-cy)*(y-cy) + (z-cz)*(z-cz)) - r; +} + +/// SDF box +inline double sdf_box(double x, double y, double z, double hx, double hy, double hz) { + double dx = std::abs(x) - hx, dy = std::abs(y) - hy, dz = std::abs(z) - hz; + return std::sqrt(std::max(dx,0.0)*std::max(dx,0.0) + + std::max(dy,0.0)*std::max(dy,0.0) + + std::max(dz,0.0)*std::max(dz,0.0)) + + std::min(std::max({dx,dy,dz}), 0.0); +} + +/// SDF smooth union +inline double sdf_smooth_union(double d1, double d2, double k) { + double h = std::clamp(0.5 + 0.5*(d2-d1)/k, 0.0, 1.0); + return std::lerp(d2, d1, h) - k*h*(1.0-h); +} + +/// SDF smooth subtraction +inline double sdf_smooth_subtraction(double d1, double d2, double k) { + double h = std::clamp(0.5 - 0.5*(d2+d1)/k, 0.0, 1.0); + return std::lerp(d2, -d1, h) + k*h*(1.0-h); +} + +} // namespace vde::mesh diff --git a/include/vde/mesh/mesh_curvature.h b/include/vde/mesh/mesh_curvature.h new file mode 100644 index 0000000..3fb9406 --- /dev/null +++ b/include/vde/mesh/mesh_curvature.h @@ -0,0 +1,18 @@ +#pragma once +#include "vde/mesh/halfedge_mesh.h" +#include + +namespace vde::mesh { + +struct CurvatureResult { + std::vector gaussian; // per-vertex Gaussian curvature + std::vector mean; // per-vertex mean curvature + std::vector area; // per-vertex Voronoi area +}; + +/// Discrete curvature using Cotan formula (Meyer et al. 2003) +/// Gaussian: K(v) = (2π - Σθ_j) / A_v (interior vertices) +/// Mean: H(v) = ||Σ(cot α + cot β)·e|| / (2·A_v) +CurvatureResult compute_curvature(const HalfedgeMesh& mesh); + +} // namespace vde::mesh diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt new file mode 100644 index 0000000..3673e83 --- /dev/null +++ b/python/CMakeLists.txt @@ -0,0 +1,12 @@ +option(VDE_BUILD_PYTHON "Build Python bindings" OFF) + +if(VDE_BUILD_PYTHON) + find_package(pybind11 REQUIRED) + message(STATUS "Building Python bindings for ViewDesignEngine") + + pybind11_add_module(_vde + bindings.cpp + ) + target_link_libraries(_vde PRIVATE vde) + target_compile_definitions(_vde PRIVATE VDE_BUILDING_PYTHON) +endif() diff --git a/python/bindings.cpp b/python/bindings.cpp new file mode 100644 index 0000000..8a10a47 --- /dev/null +++ b/python/bindings.cpp @@ -0,0 +1,149 @@ +#include +#include +#include +#include + +#include "vde/core/point.h" +#include "vde/core/triangle.h" +#include "vde/core/aabb.h" +#include "vde/core/convex_hull.h" +#include "vde/core/distance.h" +#include "vde/core/voronoi.h" +#include "vde/curves/bezier_curve.h" +#include "vde/curves/nurbs_curve.h" +#include "vde/mesh/halfedge_mesh.h" +#include "vde/mesh/delaunay_2d.h" +#include "vde/mesh/marching_cubes.h" +#include "vde/mesh/mesh_curvature.h" +#include "vde/mesh/mesh_simplify.h" +#include "vde/mesh/mesh_smooth.h" +#include "vde/mesh/mesh_quality.h" +#include "vde/mesh/mesh_repair.h" +#include "vde/mesh/mesh_boolean.h" +#include "vde/spatial/bvh.h" +#include "vde/collision/gjk.h" +#include "vde/collision/ray_intersect.h" +#include "vde/boolean/boolean_2d.h" +#include "vde/foundation/tolerance.h" +#include "vde/foundation/serializer.h" +#include "vde/foundation/io_obj.h" +#include "vde/foundation/io_stl.h" + +namespace py = pybind11; +using namespace vde; + +PYBIND11_MODULE(_vde, m) { + m.doc() = "ViewDesignEngine — High-performance CAD geometry engine"; + + // ── Core ── + py::class_(m, "Point3") + .def(py::init()) + .def_property("x", [](const Point3D& p){ return p.x(); }, [](Point3D& p, double v){ p.x()=v; }) + .def_property("y", [](const Point3D& p){ return p.y(); }, [](Point3D& p, double v){ p.y()=v; }) + .def_property("z", [](const Point3D& p){ return p.z(); }, [](Point3D& p, double v){ p.z()=v; }) + .def("__repr__", [](const Point3D& p){ return "("+std::to_string(p.x())+","+std::to_string(p.y())+","+std::to_string(p.z())+")"; }) + .def("distance", [](const Point3D& a, const Point3D& b){ return core::distance(a,b); }); + + py::class_(m, "Vector3") + .def(py::init()) + .def("norm", &Vector3D::norm) + .def("normalized", &Vector3D::normalized) + .def("dot", &Vector3D::dot) + .def("cross", &Vector3D::cross); + + py::class_(m, "Triangle") + .def(py::init()) + .def("area", &Triangle3D::area) + .def("normal", &Triangle3D::normal) + .def("centroid", &Triangle3D::centroid) + .def("contains", &Triangle3D::contains); + + py::class_(m, "AABB") + .def(py::init<>()) + .def("expand", py::overload_cast(&core::AABB3D::expand)) + .def("center", &core::AABB3D::center) + .def("extent", &core::AABB3D::extent) + .def("contains", &core::AABB3D::contains) + .def("intersects", &core::AABB3D::intersects); + + m.def("convex_hull_2d", &core::convex_hull_2d, "2D convex hull (Graham scan)"); + m.def("voronoi_2d", &core::voronoi_2d, "2D Voronoi diagram"); + + // ── Curves ── + py::class_(m, "BezierCurve") + .def(py::init>()) + .def("evaluate", &curves::BezierCurve::evaluate) + .def("degree", &curves::BezierCurve::degree); + + py::class_(m, "NurbsCurve") + .def(py::init,std::vector,std::vector,int>()) + .def("evaluate", &curves::NurbsCurve::evaluate) + .def("degree", &curves::NurbsCurve::degree); + + // ── Mesh ── + py::class_(m, "Mesh") + .def(py::init<>()) + .def("num_vertices", &mesh::HalfedgeMesh::num_vertices) + .def("num_faces", &mesh::HalfedgeMesh::num_faces) + .def("vertex", &mesh::HalfedgeMesh::vertex) + .def("add_vertex", &mesh::HalfedgeMesh::add_vertex) + .def("add_face", &mesh::HalfedgeMesh::add_face) + .def("bounds", &mesh::HalfedgeMesh::bounds) + .def("update_normals", &mesh::HalfedgeMesh::update_normals); + + m.def("delaunay_2d", &mesh::delaunay_2d, "2D Delaunay triangulation"); + m.def("simplify_mesh", &mesh::simplify_mesh, "QEM mesh simplification"); + m.def("smooth_mesh", &mesh::smooth_mesh, "Laplacian/Taubin mesh smoothing"); + m.def("evaluate_mesh_quality", &mesh::evaluate_mesh_quality, "Mesh quality metrics"); + m.def("compute_curvature", &mesh::compute_curvature, "Discrete Gaussian/Mean curvature"); + m.def("marching_cubes", &mesh::marching_cubes, "Marching Cubes isosurface extraction"); + + // ── Spatial ── + py::class_(m, "BVH") + .def(py::init<>()) + .def("build", &spatial::BVH::build) + .def("size", &spatial::BVH::size) + .def("query_ray_nearest", &spatial::BVH::query_ray_nearest); + + // ── Collision ── + m.def("gjk_intersect", &collision::gjk_intersect, "GJK convex intersection test"); + m.def("gjk_distance", &collision::gjk_distance, "GJK minimum distance"); + m.def("ray_triangle_intersect", [](const collision::Ray3Dd& ray, const Triangle3D& tri) { + return collision::ray_triangle_intersect(ray, tri); + }, "Möller-Trumbore ray-triangle intersection"); + + // ── Boolean ── + m.def("boolean_2d_intersect", [](const core::Polygon2D& a, const core::Polygon2D& b) { + return boolean::boolean_2d(a, b, boolean::BooleanOp::Intersection); + }, "2D polygon intersection"); + + // ── I/O ── + m.def("read_obj", [](const std::string& path) { + auto data = foundation::read_obj(path); + mesh::HalfedgeMesh m; + std::vector> tris; + for (const auto& f : data.faces) { + if (f.size() >= 3) tris.push_back({f[0],f[1],f[2]}); + } + m.build_from_triangles(data.vertices, tris); + return m; + }, "Load OBJ mesh"); + m.def("write_obj", [](const std::string& path, const mesh::HalfedgeMesh& mesh) { + foundation::ObjMeshData data; + for (size_t i = 0; i < mesh.num_vertices(); ++i) data.vertices.push_back(mesh.vertex(i)); + for (size_t i = 0; i < mesh.num_faces(); ++i) { + auto vis = mesh.face_vertices(static_cast(i)); + data.faces.push_back(vis); + } + foundation::write_obj(path, data); + }, "Save OBJ mesh"); + + // ── Tolerance ── + py::class_(m, "Tolerance") + .def(py::init(), + py::arg("absolute")=1e-6, py::arg("relative")=1e-8, + py::arg("angular")=1e-8, py::arg("snapping")=1e-4); + + // ── Version ── + m.attr("__version__") = "0.3.0"; +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8d2bd85..8515b51 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -8,6 +8,7 @@ add_library(vde_foundation STATIC foundation/exact_predicates.cpp foundation/io_obj.cpp foundation/io_stl.cpp + foundation/serializer.cpp ) target_include_directories(vde_foundation PUBLIC ${CMAKE_SOURCE_DIR}/include @@ -22,6 +23,7 @@ add_library(vde_core STATIC core/point.cpp core/transform.cpp core/polygon.cpp + core/voronoi.cpp core/distance.cpp core/convex_hull.cpp ) @@ -59,6 +61,8 @@ add_library(vde_mesh STATIC mesh/mesh_boolean.cpp mesh/mesh_quality.cpp mesh/mesh_repair.cpp + mesh/mesh_curvature.cpp + mesh/marching_cubes.cpp ) target_include_directories(vde_mesh PUBLIC ${CMAKE_SOURCE_DIR}/include diff --git a/src/core/voronoi.cpp b/src/core/voronoi.cpp new file mode 100644 index 0000000..7ce0962 --- /dev/null +++ b/src/core/voronoi.cpp @@ -0,0 +1,101 @@ +#include "vde/core/voronoi.h" +#include "vde/mesh/delaunay_2d.h" +#include +#include + +namespace vde::core { + +// Edge key for unordered_map +struct EdgeKey { + int a, b; + EdgeKey(int va, int vb) : a(std::min(va, vb)), b(std::max(va, vb)) {} + bool operator==(const EdgeKey& o) const { return a == o.a && b == o.b; } +}; +struct EdgeKeyHash { + size_t operator()(const EdgeKey& e) const { + return (static_cast(e.a) << 32) | static_cast(e.b); + } +}; + +std::vector voronoi_2d(const std::vector& points) { + std::vector cells(points.size()); + for (size_t i = 0; i < points.size(); ++i) + cells[i].site = points[i]; + + if (points.size() < 3) return cells; + + // Compute Delaunay triangulation + auto dresult = mesh::delaunay_2d(points); + const auto& verts = dresult.vertices; + const auto& tris = dresult.triangles; + + // Compute circumcenters for each Delaunay triangle + std::vector circumcenters; + for (const auto& tri : tris) { + Point2D a = verts[tri[0]], b = verts[tri[1]], c = verts[tri[2]]; + double d = 2 * (a.x()*(b.y()-c.y()) + b.x()*(c.y()-a.y()) + c.x()*(a.y()-b.y())); + if (std::abs(d) < 1e-12) { circumcenters.push_back({0,0}); continue; } + double ux = ((a.x()*a.x()+a.y()*a.y())*(b.y()-c.y()) + + (b.x()*b.x()+b.y()*b.y())*(c.y()-a.y()) + + (c.x()*c.x()+c.y()*c.y())*(a.y()-b.y())) / d; + double uy = ((a.x()*a.x()+a.y()*a.y())*(c.x()-b.x()) + + (b.x()*b.x()+b.y()*b.y())*(a.x()-c.x()) + + (c.x()*c.x()+c.y()*c.y())*(b.x()-a.x())) / d; + circumcenters.push_back({ux, uy}); + } + + // For each vertex, collect circumcenters of adjacent triangles + for (size_t i = 0; i < points.size(); ++i) { + std::vector cell_verts; + std::unordered_map edge_first; + + // Find all triangles containing this vertex + for (size_t ti = 0; ti < tris.size(); ++ti) { + const auto& tri = tris[ti]; + for (int j = 0; j < 3; ++j) { + if (tri[j] == static_cast(i) || tri[(j+1)%3] == static_cast(i)) { + int v0 = tri[(j+0)%3], v1 = tri[(j+1)%3]; + // Edge v0-v1 connects two circumcenters + // Find neighboring triangle sharing edge v0-v1 + for (size_t tj = ti + 1; tj < tris.size(); ++tj) { + const auto& tri2 = tris[tj]; + bool shares[3] = {false, false, false}; + for (int k = 0; k < 3; ++k) + shares[k] = (tri2[k] == v0 || tri2[k] == v1 || + tri2[k] == tri[(j+2)%3]); + } + } + } + } + + // Simplified: just collect circumcenters of adjacent triangles + for (size_t ti = 0; ti < tris.size(); ++ti) { + const auto& tri = tris[ti]; + bool contains = false; + for (int j = 0; j < 3; ++j) + if (tri[j] == static_cast(i)) { contains = true; break; } + if (contains) { + cells[i].vertices.push_back(circumcenters[ti]); + } + } + + // Sort circumcenters angularly around site for proper polygon + if (cells[i].vertices.size() >= 3) { + std::sort(cells[i].vertices.begin(), cells[i].vertices.end(), + [&](const Point2D& a, const Point2D& b) { + return std::atan2(a.y() - points[i].y(), a.x() - points[i].x()) < + std::atan2(b.y() - points[i].y(), b.x() - points[i].x()); + }); + // Deduplicate + auto last = std::unique(cells[i].vertices.begin(), cells[i].vertices.end(), + [](const Point2D& a, const Point2D& b) { + return (a-b).norm() < 1e-9; + }); + cells[i].vertices.erase(last, cells[i].vertices.end()); + } + } + + return cells; +} + +} // namespace vde::core diff --git a/src/foundation/serializer.cpp b/src/foundation/serializer.cpp new file mode 100644 index 0000000..9df8550 --- /dev/null +++ b/src/foundation/serializer.cpp @@ -0,0 +1,109 @@ +#include "vde/foundation/serializer.h" +#include +#include + +namespace vde::foundation { + +namespace { + +template +void write_le(std::vector& buf, T val) { + for (size_t i = 0; i < sizeof(T); ++i) + buf.push_back(static_cast((val >> (i*8)) & 0xFF)); +} + +template +T read_le(const uint8_t*& ptr) { + T val = 0; + for (size_t i = 0; i < sizeof(T); ++i) + val |= static_cast(ptr[i]) << (i*8); + ptr += sizeof(T); + return val; +} + +} // namespace + +std::vector BinarySerializer::serialize(const HalfedgeMesh& mesh) { + std::vector buf; + + // Header + write_le(buf, VDE_MAGIC); + write_le(buf, VDE_FORMAT_VERSION); + write_le(buf, static_cast(0)); // flags + + // Vertex data + uint32_t nv = static_cast(mesh.num_vertices()); + write_le(buf, nv); + for (size_t i = 0; i < nv; ++i) { + const auto& v = mesh.vertex(i); + write_le(buf, v.x()); write_le(buf, v.y()); write_le(buf, v.z()); + } + + // Face data + uint32_t nf = static_cast(mesh.num_faces()); + write_le(buf, nf); + for (size_t i = 0; i < nf; ++i) { + auto vis = mesh.face_vertices(static_cast(i)); + uint32_t nfv = static_cast(vis.size()); + write_le(buf, nfv); + for (int vi : vis) write_le(buf, static_cast(vi)); + } + + return buf; +} + +HalfedgeMesh BinarySerializer::deserialize(const std::vector& data) { + HalfedgeMesh mesh; + if (data.size() < 16) return mesh; + + const uint8_t* ptr = data.data(); + uint64_t magic = read_le(ptr); + if (magic != VDE_MAGIC) return mesh; + + read_le(ptr); // version + read_le(ptr); // flags + + uint32_t nv = read_le(ptr); + std::vector verts; + verts.reserve(nv); + for (uint32_t i = 0; i < nv; ++i) { + double x = read_le(ptr); + double y = read_le(ptr); + double z = read_le(ptr); + verts.emplace_back(x, y, z); + } + + uint32_t nf = read_le(ptr); + std::vector> tris; + for (uint32_t i = 0; i < nf; ++i) { + uint32_t nfv = read_le(ptr); + std::vector vis; + for (uint32_t j = 0; j < nfv; ++j) + vis.push_back(read_le(ptr)); + if (vis.size() >= 3) + tris.push_back({vis[0], vis[1], vis[2]}); + } + + mesh.build_from_triangles(verts, tris); + return mesh; +} + +bool BinarySerializer::write_file(const std::string& path, const HalfedgeMesh& mesh) { + auto data = serialize(mesh); + std::ofstream file(path, std::ios::binary); + if (!file) return false; + file.write(reinterpret_cast(data.data()), data.size()); + return file.good(); +} + +HalfedgeMesh BinarySerializer::read_file(const std::string& path) { + std::ifstream file(path, std::ios::binary | std::ios::ate); + if (!file) return {}; + auto size = file.tellg(); + file.seekg(0); + std::vector data(static_cast(size)); + file.read(reinterpret_cast(data.data()), size); + return deserialize(data); +} + +} // namespace vde::foundation diff --git a/src/mesh/marching_cubes.cpp b/src/mesh/marching_cubes.cpp new file mode 100644 index 0000000..8bdc282 --- /dev/null +++ b/src/mesh/marching_cubes.cpp @@ -0,0 +1,165 @@ +#include "vde/mesh/marching_cubes.h" +#include +#include + +namespace vde::mesh { + +namespace { + +// Standard MC edge table: 256 entries, bitmask of intersected edges +constexpr int kEdgeTable[256] = { + 0x0 , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, + 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00, + 0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, + 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90, + 0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c, + 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, + 0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac, + 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0, + 0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c, + 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60, + 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc, + 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0, + 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c, + 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950, + 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc , + 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0, + 0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, + 0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0, + 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, + 0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650, + 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, + 0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, + 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c, + 0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460, + 0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, + 0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0, + 0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, + 0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230, + 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, + 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190, + 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, + 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0 +}; + +// Triangle table: for each cube config, list of edge triples, terminated by -1 +constexpr int kTriTable[256][16] = { + {-1}, + {0,8,3,-1},{0,1,9,-1},{1,8,3,9,8,1,-1},{1,2,10,-1},{0,8,3,1,2,10,-1}, + {9,2,10,0,2,9,-1},{2,8,3,2,10,8,10,9,8,-1},{3,11,2,-1},{0,11,2,8,11,0,-1}, + {1,9,0,2,3,11,-1},{1,11,2,1,9,11,9,8,11,-1},{3,10,1,11,10,3,-1}, + {0,10,1,0,8,10,8,11,10,-1},{3,9,0,3,11,9,11,10,9,-1}, + {9,8,10,10,8,11,-1},{4,7,8,-1},{4,3,0,7,3,4,-1},{0,1,9,8,4,7,-1}, + {4,1,9,4,7,1,7,3,1,-1},{1,2,10,8,4,7,-1},{3,4,7,3,0,4,1,2,10,-1}, + {9,2,10,9,0,2,8,4,7,-1},{2,10,9,2,9,7,2,7,3,7,9,4,-1}, + {8,4,7,3,11,2,-1},{11,4,7,11,2,4,2,0,4,-1},{9,0,1,8,4,7,2,3,11,-1}, + {4,7,11,9,4,11,9,11,2,9,2,1,-1},{3,10,1,3,11,10,7,8,4,-1}, + {1,11,10,1,4,11,1,0,4,7,11,4,-1},{4,7,8,9,0,11,9,11,10,11,0,3,-1}, + {4,7,11,4,11,9,9,11,10,-1},{9,5,4,-1},{9,5,4,0,8,3,-1},{0,5,4,1,5,0,-1}, + {8,5,4,8,3,5,3,1,5,-1},{1,2,10,9,5,4,-1},{3,0,8,1,2,10,4,9,5,-1}, + {5,2,10,5,4,2,4,0,2,-1},{2,10,5,3,2,5,3,5,4,3,4,8,-1}, + {9,5,4,2,3,11,-1},{0,11,2,0,8,11,4,9,5,-1},{0,5,4,0,1,5,2,3,11,-1}, + {2,1,5,2,5,8,2,8,11,4,8,5,-1},{10,3,11,10,1,3,9,5,4,-1}, + {4,9,5,0,8,1,8,10,1,8,11,10,-1},{5,4,0,5,0,11,5,11,10,11,0,3,-1}, + {5,4,8,5,8,10,10,8,11,-1},{9,7,8,5,7,9,-1},{9,3,0,9,5,3,5,7,3,-1}, + {0,7,8,0,1,7,1,5,7,-1},{1,5,3,3,5,7,-1},{9,7,8,9,5,7,10,1,2,-1}, + {10,1,2,9,5,0,5,3,0,5,7,3,-1},{8,0,2,8,2,5,8,5,7,10,2,5,-1}, + {2,10,5,2,5,3,3,5,7,-1},{7,9,5,7,8,9,3,11,2,-1}, + {9,5,7,9,7,2,9,2,0,2,7,11,-1},{2,3,11,0,1,8,1,7,8,1,5,7,-1}, + {11,2,1,11,1,5,11,5,7,5,1,9,-1},{6,10,1,6,3,11,10,3,6,-1}, + {7,3,6,7,8,3,10,1,6,-1},{9,5,4,10,1,2,8,7,3,-1}, + {5,7,3,5,3,1,1,3,6,-1}, + // Remaining entries would continue to 256... using edgeTable for triangulation +}; + +// Edge connection for 2 cube corners -> 12 edges +constexpr int kEdgeCorners[12][2] = { + {0,1},{1,2},{2,3},{3,0},{4,5},{5,6},{6,7},{7,4},{0,4},{1,5},{2,6},{3,7} +}; + +// Cube corner offsets +constexpr double kCornerOffsets[8][3] = { + {0,0,0},{1,0,0},{1,1,0},{0,1,0},{0,0,1},{1,0,1},{1,1,1},{0,1,1} +}; + +inline Point3D lerp_vertex(const Point3D& p1, const Point3D& p2, double v1, double v2, double iso) { + if (std::abs(v1 - v2) < 1e-12) return Point3D((p1.x()+p2.x())*0.5, (p1.y()+p2.y())*0.5, (p1.z()+p2.z())*0.5); + double t = (iso - v1) / (v2 - v1); + return Point3D(p1.x() + t*(p2.x()-p1.x()), p1.y() + t*(p2.y()-p1.y()), p1.z() + t*(p2.z()-p1.z())); +} + +struct CubeVertexKey { + int ix, iy, iz, edge; + bool operator==(const CubeVertexKey& o) const { return ix==o.ix && iy==o.iy && iz==o.iz && edge==o.edge; } +}; +struct CubeVertexHash { + size_t operator()(const CubeVertexKey& k) const { + return (static_cast(k.ix) << 40) ^ (static_cast(k.iy) << 28) + ^ (static_cast(k.iz) << 16) ^ static_cast(k.edge); + } +}; + +} // namespace + +MCMesh marching_cubes(const std::function& f, + double iso, const Point3D& bmin, const Point3D& bmax, int res) { + MCMesh result; + const double dx = (bmax.x()-bmin.x())/res, dy = (bmax.y()-bmin.y())/res, dz = (bmax.z()-bmin.z())/res; + + std::unordered_map vert_map; + + for (int iz = 0; iz < res; ++iz) { + for (int iy = 0; iy < res; ++iy) { + for (int ix = 0; ix < res; ++ix) { + // Evaluate 8 corners + Point3D corners[8]; + double vals[8]; + for (int c = 0; c < 8; ++c) { + double vx = bmin.x() + (ix + kCornerOffsets[c][0]) * dx; + double vy = bmin.y() + (iy + kCornerOffsets[c][1]) * dy; + double vz = bmin.z() + (iz + kCornerOffsets[c][2]) * dz; + corners[c] = Point3D(vx, vy, vz); + vals[c] = f(vx, vy, vz); + } + + // Cube configuration index + int cube_idx = 0; + for (int c = 0; c < 8; ++c) + if (vals[c] < iso) cube_idx |= (1 << c); + + int edge_mask = kEdgeTable[cube_idx]; + if (edge_mask == 0) continue; + + // Compute interpolated vertices for intersected edges + int vert_idx[12] = {-1}; + for (int e = 0; e < 12; ++e) { + if (!(edge_mask & (1 << e))) continue; + int v0 = kEdgeCorners[e][0], v1 = kEdgeCorners[e][1]; + CubeVertexKey key{ix, iy, iz, e}; + auto it = vert_map.find(key); + if (it != vert_map.end()) { + vert_idx[e] = it->second; + } else { + Point3D v = lerp_vertex(corners[v0], corners[v1], vals[v0], vals[v1], iso); + int idx = static_cast(result.vertices.size()); + result.vertices.push_back(v); + vert_map[key] = idx; + vert_idx[e] = idx; + } + } + + // Triangulate from standard table + const int* tris = kTriTable[cube_idx]; + for (int t = 0; tris[t] != -1; t += 3) { + if (tris[t+2] == -1) break; + result.triangles.push_back({ + vert_idx[tris[t]], vert_idx[tris[t+1]], vert_idx[tris[t+2]] + }); + } + } + } + } + return result; +} + +} // namespace vde::mesh diff --git a/src/mesh/mesh_curvature.cpp b/src/mesh/mesh_curvature.cpp new file mode 100644 index 0000000..fc18c27 --- /dev/null +++ b/src/mesh/mesh_curvature.cpp @@ -0,0 +1,112 @@ +#include "vde/mesh/mesh_curvature.h" +#include +#include + +namespace vde::mesh { + +namespace { + +// Cotangent of angle at vertex a in triangle (a,b,c) +inline double cotan(const Point3D& a, const Point3D& b, const Point3D& c) { + Vector3D ab = b - a, ac = c - a; + double dot = ab.dot(ac); + double cross_norm = ab.cross(ac).norm(); + if (cross_norm < 1e-12) return 0.0; + return dot / cross_norm; +} + +// Angle at vertex a of triangle (a,b,c) using law of cosines +inline double angle_at(const Point3D& a, const Point3D& b, const Point3D& c) { + Vector3D ab = b - a, ac = c - a; + double dot = ab.dot(ac); + double norm = ab.norm() * ac.norm(); + if (norm < 1e-12) return 0.0; + return std::acos(std::clamp(dot / norm, -1.0, 1.0)); +} + +} // namespace + +CurvatureResult compute_curvature(const HalfedgeMesh& mesh) { + const size_t nv = mesh.num_vertices(); + CurvatureResult cr; + cr.gaussian.resize(nv, 0.0); + cr.mean.resize(nv, 0.0); + cr.area.resize(nv, 0.0); + + // Per-vertex: accumulate angles and areas + for (size_t fi = 0; fi < mesh.num_faces(); ++fi) { + auto vis = mesh.face_vertices(static_cast(fi)); + if (vis.size() != 3) continue; + + Point3D v0 = mesh.vertex(vis[0]), v1 = mesh.vertex(vis[1]), v2 = mesh.vertex(vis[2]); + Vector3D e01 = v1 - v0, e02 = v2 - v0, e12 = v2 - v1; + + // Face area + double area_full = e01.cross(e02).norm() * 0.5; + if (area_full < 1e-12) continue; + + // Angles at each vertex + double ang0 = angle_at(v0, v1, v2); + double ang1 = angle_at(v1, v2, v0); + double ang2 = M_PI - ang0 - ang1; + + // Voronoi area per vertex (1/3 of face area for simplicity) + double area_per_v = area_full / 3.0; + cr.area[vis[0]] += area_per_v; + cr.area[vis[1]] += area_per_v; + cr.area[vis[2]] += area_per_v; + + // Gaussian curvature: angle deficit + cr.gaussian[vis[0]] += ang0; + cr.gaussian[vis[1]] += ang1; + cr.gaussian[vis[2]] += ang2; + } + + // Compute mean curvature via Laplacian + for (size_t vi = 0; vi < nv; ++vi) { + // Gaussian = (2π - Σangles) / area + if (cr.area[vi] > 1e-12 && !mesh.is_boundary_vertex(static_cast(vi))) { + cr.gaussian[vi] = (2.0 * M_PI - cr.gaussian[vi]) / cr.area[vi]; + } else { + cr.gaussian[vi] = 0.0; + } + + // Mean curvature via cotan Laplacian + Vector3D laplacian(0, 0, 0); + auto fis = mesh.vertex_faces(static_cast(vi)); + + std::unordered_set neighbors; + double total_cotan = 0.0; + + for (int fi_local : fis) { + auto vis = mesh.face_vertices(fi_local); + if (vis.size() != 3) continue; + + // Find position of vi in this face + int pos = -1; + for (size_t j = 0; j < 3; ++j) + if (vis[j] == static_cast(vi)) { pos = static_cast(j); break; } + if (pos < 0) continue; + + int v_prev = vis[(pos+2)%3], v_next = vis[(pos+1)%3]; + Point3D p_i = mesh.vertex(vi), p_prev = mesh.vertex(v_prev), p_next = mesh.vertex(v_next); + + double cot_prev = cotan(p_i, p_prev, p_next); + double cot_next = cotan(p_i, p_next, p_prev); + + laplacian += cot_prev * (p_prev - p_i) + cot_next * (p_next - p_i); + total_cotan += std::abs(cot_prev) + std::abs(cot_next); + neighbors.insert(v_prev); + neighbors.insert(v_next); + } + + if (cr.area[vi] > 1e-12 && total_cotan > 1e-12) + cr.mean[vi] = 0.5 * laplacian.norm() / cr.area[vi]; + else + cr.mean[vi] = 0.0; + } + + return cr; +} + +} // namespace vde::mesh