feat: v0.5.0 P1 + 差异化功能 — Voronoi, Marching Cubes, 曲率, 序列化, Python绑定
新增功能： - core: 2D Voronoi 图（Delaunay 对偶图） - mesh: Marching Cubes（完整256项三角表 + SDF 辅助函数） - mesh: 离散曲率（Gaussian/Mean, Cotan 公式） - foundation: 二进制序列化（版本化格式，Little-Endian） - python: pybind11 绑定（30+ API 暴露） - marching_cubes.h: SDF 基本体（球/盒）+ 光滑布尔
This commit is contained in:
ViewDesignEngine
@@ -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
|
||||
|
||||
@@ -0,0 +1,101 @@
|
||||
#include "vde/core/voronoi.h"
|
||||
#include "vde/mesh/delaunay_2d.h"
|
||||
#include <unordered_map>
|
||||
#include <algorithm>
|
||||
|
||||
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<uint64_t>(e.a) << 32) | static_cast<uint64_t>(e.b);
|
||||
}
|
||||
};
|
||||
|
||||
std::vector<VoronoiCell> voronoi_2d(const std::vector<Point2D>& points) {
|
||||
std::vector<VoronoiCell> 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<Point2D> 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<Point2D> cell_verts;
|
||||
std::unordered_map<EdgeKey, int, EdgeKeyHash> 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<int>(i) || tri[(j+1)%3] == static_cast<int>(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<int>(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
|
||||
@@ -0,0 +1,109 @@
|
||||
#include "vde/foundation/serializer.h"
|
||||
#include <fstream>
|
||||
#include <cstring>
|
||||
|
||||
namespace vde::foundation {
|
||||
|
||||
namespace {
|
||||
|
||||
template<typename T>
|
||||
void write_le(std::vector<uint8_t>& buf, T val) {
|
||||
for (size_t i = 0; i < sizeof(T); ++i)
|
||||
buf.push_back(static_cast<uint8_t>((val >> (i*8)) & 0xFF));
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
T read_le(const uint8_t*& ptr) {
|
||||
T val = 0;
|
||||
for (size_t i = 0; i < sizeof(T); ++i)
|
||||
val |= static_cast<T>(ptr[i]) << (i*8);
|
||||
ptr += sizeof(T);
|
||||
return val;
|
||||
}
|
||||
|
||||
} // namespace
|
||||
|
||||
std::vector<uint8_t> BinarySerializer::serialize(const HalfedgeMesh& mesh) {
|
||||
std::vector<uint8_t> buf;
|
||||
|
||||
// Header
|
||||
write_le(buf, VDE_MAGIC);
|
||||
write_le(buf, VDE_FORMAT_VERSION);
|
||||
write_le(buf, static_cast<uint32_t>(0)); // flags
|
||||
|
||||
// Vertex data
|
||||
uint32_t nv = static_cast<uint32_t>(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<uint32_t>(mesh.num_faces());
|
||||
write_le(buf, nf);
|
||||
for (size_t i = 0; i < nf; ++i) {
|
||||
auto vis = mesh.face_vertices(static_cast<int>(i));
|
||||
uint32_t nfv = static_cast<uint32_t>(vis.size());
|
||||
write_le(buf, nfv);
|
||||
for (int vi : vis) write_le(buf, static_cast<int32_t>(vi));
|
||||
}
|
||||
|
||||
return buf;
|
||||
}
|
||||
|
||||
HalfedgeMesh BinarySerializer::deserialize(const std::vector<uint8_t>& data) {
|
||||
HalfedgeMesh mesh;
|
||||
if (data.size() < 16) return mesh;
|
||||
|
||||
const uint8_t* ptr = data.data();
|
||||
uint64_t magic = read_le<uint64_t>(ptr);
|
||||
if (magic != VDE_MAGIC) return mesh;
|
||||
|
||||
read_le<uint32_t>(ptr); // version
|
||||
read_le<uint32_t>(ptr); // flags
|
||||
|
||||
uint32_t nv = read_le<uint32_t>(ptr);
|
||||
std::vector<Point3D> verts;
|
||||
verts.reserve(nv);
|
||||
for (uint32_t i = 0; i < nv; ++i) {
|
||||
double x = read_le<double>(ptr);
|
||||
double y = read_le<double>(ptr);
|
||||
double z = read_le<double>(ptr);
|
||||
verts.emplace_back(x, y, z);
|
||||
}
|
||||
|
||||
uint32_t nf = read_le<uint32_t>(ptr);
|
||||
std::vector<std::array<int, 3>> tris;
|
||||
for (uint32_t i = 0; i < nf; ++i) {
|
||||
uint32_t nfv = read_le<uint32_t>(ptr);
|
||||
std::vector<int> vis;
|
||||
for (uint32_t j = 0; j < nfv; ++j)
|
||||
vis.push_back(read_le<int32_t>(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<const char*>(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<uint8_t> data(static_cast<size_t>(size));
|
||||
file.read(reinterpret_cast<char*>(data.data()), size);
|
||||
return deserialize(data);
|
||||
}
|
||||
|
||||
} // namespace vde::foundation
|
||||
@@ -0,0 +1,165 @@
|
||||
#include "vde/mesh/marching_cubes.h"
|
||||
#include <unordered_map>
|
||||
#include <cmath>
|
||||
|
||||
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<uint64_t>(k.ix) << 40) ^ (static_cast<uint64_t>(k.iy) << 28)
|
||||
^ (static_cast<uint64_t>(k.iz) << 16) ^ static_cast<uint64_t>(k.edge);
|
||||
}
|
||||
};
|
||||
|
||||
} // namespace
|
||||
|
||||
MCMesh marching_cubes(const std::function<double(double,double,double)>& 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<CubeVertexKey, int, CubeVertexHash> 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<int>(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
|
||||
@@ -0,0 +1,112 @@
|
||||
#include "vde/mesh/mesh_curvature.h"
|
||||
#include <unordered_set>
|
||||
#include <cmath>
|
||||
|
||||
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<int>(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<int>(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<int>(vi));
|
||||
|
||||
std::unordered_set<int> 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<int>(vi)) { pos = static_cast<int>(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
|
||||
Reference in New Issue
Block a user