11#include <unordered_set>
14#include <igl/read_triangle_mesh.h>
16#include <geogram/mesh/mesh_io.h>
17#include <geogram/mesh/mesh_geometry.h>
18#include <geogram/basic/geometry.h>
19#include <geogram/mesh/mesh_preprocessing.h>
20#include <geogram/mesh/mesh_topology.h>
21#include <geogram/mesh/mesh_geometry.h>
22#include <geogram/mesh/mesh_repair.h>
23#include <geogram/mesh/mesh_AABB.h>
24#include <geogram/voronoi/CVT.h>
25#include <geogram/basic/logger.h>
33 if (M.vertices.dimension() == 2)
36 assert(M.vertices.dimension() == 3);
37 GEO::vec3 min_corner, max_corner;
38 GEO::get_bbox(M, &min_corner[0], &max_corner[0]);
39 const double diff = (max_corner[2] - min_corner[2]);
41 return fabs(diff) < tol;
48 for (index_t d = 0; d < std::min(3u, (index_t)M.vertices.dimension()); ++d)
50 if (M.vertices.double_precision())
52 p[d] = M.vertices.point_ptr(v)[d];
56 p[d] = M.vertices.single_precision_point_ptr(v)[d];
68 for (index_t lv = 0; lv < M.facets.nb_vertices(f); ++lv)
72 return p / M.facets.nb_vertices(f);
80 auto v = M.vertices.create_vertex();
81 for (index_t d = 0; d < std::min(3u, (index_t)M.vertices.dimension()); ++d)
83 if (M.vertices.double_precision())
85 M.vertices.point_ptr(v)[d] = p[d];
89 M.vertices.single_precision_point_ptr(v)[d] = (float)p[d];
101 std::vector<ElementType> old_tags = element_tags;
103 element_tags.resize(M.facets.nb());
106 GEO::Attribute<bool> is_boundary_vertex(M.vertices.attributes(),
"boundary_vertex");
107 std::vector<bool> is_interface_vertex(M.vertices.nb(),
false);
109 for (index_t f = 0; f < M.facets.nb(); ++f)
111 if (M.facets.nb_vertices(f) != 4 || (!old_tags.empty() && old_tags[f] == ElementType::INTERIOR_POLYTOPE))
114 for (index_t lv = 0; lv < M.facets.nb_vertices(f); ++lv)
116 is_interface_vertex[M.facets.vertex(f, lv)] =
true;
126 std::vector<int> degree(M.vertices.nb(), 0);
127 std::vector<bool> is_regular_vertex(M.vertices.nb());
128 for (index_t f = 0; f < M.facets.nb(); ++f)
130 if (M.facets.nb_vertices(f) == 4)
133 for (index_t lv = 0; lv < M.facets.nb_vertices(f); ++lv)
135 index_t v = M.facets.vertex(f, lv);
140 for (index_t v = 0; v < M.vertices.nb(); ++v)
143 if (is_boundary_vertex[v] || is_interface_vertex[v])
145 is_regular_vertex[v] = (degree[v] <= 2);
149 is_regular_vertex[v] = (degree[v] == 4);
154 for (index_t f = 0; f < M.facets.nb(); ++f)
156 assert(M.facets.nb_vertices(f) > 2);
157 if (!old_tags.empty() && old_tags[f] == ElementType::INTERIOR_POLYTOPE)
160 if (M.facets.nb_vertices(f) == 4)
165 bool is_boundary_facet =
false;
166 bool is_interface_facet =
false;
167 for (index_t lv = 0; lv < M.facets.nb_vertices(f); ++lv)
169 if (is_boundary_vertex[M.facets.vertex(f, lv)])
171 is_boundary_facet =
true;
173 if (is_interface_vertex[M.facets.vertex(f, lv)])
175 is_interface_facet =
true;
180 if (is_boundary_facet || is_interface_facet)
184 bool is_singular =
false;
185 for (index_t lv = 0; lv < M.facets.nb_vertices(f); ++lv)
187 index_t v = M.facets.vertex(f, lv);
188 if (is_boundary_vertex[v] || is_interface_vertex[v])
190 if (!is_regular_vertex[v])
198 if (!is_regular_vertex[v])
200 element_tags[f] = ElementType::UNDEFINED;
206 if (is_interface_facet)
208 element_tags[f] = ElementType::INTERFACE_CUBE;
210 else if (is_singular)
212 element_tags[f] = ElementType::SIMPLE_SINGULAR_BOUNDARY_CUBE;
216 element_tags[f] = ElementType::REGULAR_BOUNDARY_CUBE;
222 int nb_singulars = 0;
223 for (index_t lv = 0; lv < M.facets.nb_vertices(f); ++lv)
225 if (!is_regular_vertex[M.facets.vertex(f, lv)])
231 if (nb_singulars == 0)
233 element_tags[f] = ElementType::REGULAR_INTERIOR_CUBE;
235 else if (nb_singulars == 1)
237 element_tags[f] = ElementType::SIMPLE_SINGULAR_INTERIOR_CUBE;
241 element_tags[f] = ElementType::MULTI_SINGULAR_INTERIOR_CUBE;
251 GEO::Attribute<bool> boundary_vertices(M.vertices.attributes(),
"boundary_vertex");
252 for (index_t lv = 0; lv < M.facets.nb_vertices(f); ++lv)
254 if (boundary_vertices[M.facets.vertex(f, lv)])
256 tag = ElementType::BOUNDARY_POLYTOPE;
261 element_tags[f] = tag;
267 for (index_t f = 0; f < M.facets.nb(); ++f)
269 if (M.facets.nb_vertices(f) == 3)
271 element_tags[f] = ElementType::SIMPLEX;
282 double signed_area(
const GEO::Mesh &M, GEO::index_t f)
286 index_t v0 = M.facet_corners.vertex(M.facets.corners_begin(f));
287 const vec3 &p0 = Geom::mesh_vertex(M, v0);
289 M.facets.corners_begin(f) + 1;
290 c + 1 < M.facets.corners_end(f); c++)
292 index_t v1 = M.facet_corners.vertex(c);
294 index_t v2 = M.facet_corners.vertex(c + 1);
296 result += Geom::triangle_signed_area(vec2(&p0[0]), vec2(&p1[0]), vec2(&p2[0]));
306 vector<index_t> component;
307 index_t nb_components = get_connected_components(M, component);
308 vector<double> comp_signed_volume(nb_components, 0.0);
309 for (index_t f = 0; f < M.facets.nb(); ++f)
311 comp_signed_volume[component[f]] += signed_area(M, f);
313 for (index_t f = 0; f < M.facets.nb(); ++f)
315 if (comp_signed_volume[component[f]] < 0.0)
326 assert(
V.rows() == C.size());
327 int num_colors = C.maxCoeff() + 1;
328 Eigen::VectorXi count(num_colors);
330 for (
int i = 0; i < C.size(); ++i)
334 R.resize(num_colors + 1);
336 for (
int c = 0; c < num_colors; ++c)
338 R(c + 1) = R(c) + count(c);
341 Eigen::VectorXi remap(C.size());
342 for (
int i = 0; i < C.size(); ++i)
344 remap[i] = R(C(i)) + count[C(i)];
348 Eigen::MatrixXd NV(
V.rows(),
V.cols());
349 for (
int v = 0; v <
V.rows(); ++v)
351 NV.row(remap(v)) =
V.row(v);
355 for (
int f = 0; f < F.rows(); ++f)
357 for (
int lv = 0; lv < F.cols(); ++lv)
359 F(f, lv) = remap(F(f, lv));
369 void compute_unsigned_distance_field(
const GEO::Mesh &M,
370 const GEO::MeshFacetsAABB &aabb_tree,
const Eigen::MatrixXd &P, Eigen::VectorXd &D)
372 assert(P.cols() == 3);
374#pragma omp parallel for
375 for (
int i = 0; i < P.rows(); ++i)
377 GEO::vec3 pos(P(i, 0), P(i, 1), P(i, 2));
378 double sq_dist = aabb_tree.squared_distance(pos);
386 double x1,
double y1,
double x2,
double y2,
double &twice_signed_area)
388 twice_signed_area = y1 * x2 - x1 * y2;
389 if (twice_signed_area > 0)
391 else if (twice_signed_area < 0)
410 bool point_in_triangle_2d(
411 double x0,
double y0,
double x1,
double y1,
412 double x2,
double y2,
double x3,
double y3,
413 double &a,
double &b,
double &c)
421 int signa = orientation(x2, y2, x3, y3, a);
424 int signb = orientation(x3, y3, x1, y1, b);
427 int signc = orientation(x1, y1, x2, y2, c);
430 double sum = a + b + c;
431 geo_assert(sum != 0);
452 inline GEO::Sign orient_2d_inexact(GEO::vec2 p0, GEO::vec2 p1, GEO::vec2 p2)
454 double a11 = p1[0] - p0[0];
455 double a12 = p1[1] - p0[1];
457 double a21 = p2[0] - p0[0];
458 double a22 = p2[1] - p0[1];
460 double Delta = GEO::det2x2(
464 return GEO::geo_sgn(Delta);
479 template <
int X = 0,
int Y = 1,
int Z = 2>
480 int intersect_ray_z(
const GEO::Mesh &M, GEO::index_t f,
const GEO::vec3 &q,
double &
z)
484 index_t c =
M.facets.corners_begin(f);
485 const vec3 &p1 = Geom::mesh_vertex(M,
M.facet_corners.vertex(c++));
486 const vec3 &p2 = Geom::mesh_vertex(M,
M.facet_corners.vertex(c++));
487 const vec3 &p3 = Geom::mesh_vertex(M,
M.facet_corners.vertex(c));
490 if (point_in_triangle_2d(
491 q[X], q[Y], p1[X], p1[Y], p2[X], p2[Y], p3[X], p3[Y], u, v, w))
493 z = u * p1[Z] + v * p2[Z] + w * p3[Z];
494 auto sign = orient_2d_inexact(vec2(p1[X], p1[Y]), vec2(p2[X], p2[Y]), vec2(p3[X], p3[Y]));
512 void compute_sign(
const GEO::Mesh &M,
const GEO::MeshFacetsAABB &aabb_tree,
513 const Eigen::MatrixXd &P, Eigen::VectorXd &D)
515 assert(
P.cols() == 3);
516 assert(D.size() ==
P.rows());
518 GEO::vec3 min_corner, max_corner;
519 GEO::get_bbox(M, &min_corner[0], &max_corner[0]);
521#pragma omp parallel for
522 for (
int k = 0; k <
P.rows(); ++k)
524 GEO::vec3 center(
P(k, 0),
P(k, 1),
P(k, 2));
527 box.xyz_min[0] = box.xyz_max[0] = center[0];
528 box.xyz_min[1] = box.xyz_max[1] = center[1];
529 box.xyz_min[2] = min_corner[2];
530 box.xyz_max[2] = max_corner[2];
532 std::vector<std::pair<double, int>> inter;
533 auto action = [&
M, &inter, ¢er](GEO::index_t
f) {
535 if (
int s = intersect_ray_z(M, f, center,
z))
537 inter.emplace_back(
z, s);
540 aabb_tree.compute_bbox_facet_bbox_intersections(box, action);
541 std::sort(inter.begin(), inter.end());
543 std::vector<double> reduced;
544 for (
int i = 0, s = 0; i < (int)inter.size(); ++i)
546 const int ds = inter[i].second;
548 if ((s == -1 && ds < 0) || (s == 0 && ds > 0))
550 reduced.push_back(inter[i].first);
555 for (
double z : reduced)
562 if (num_before % 2 == 1)
578 M.vertices.create_vertices((
int)
V.rows());
579 for (
int i = 0; i < (int)M.vertices.nb(); ++i)
581 GEO::vec3 &p = M.vertices.point(i);
584 p[2] =
V.cols() >= 3 ?
V(i, 2) : 0;
589 M.facets.create_triangles((
int)F.rows());
591 else if (F.cols() == 4)
593 M.facets.create_quads((
int)F.rows());
597 throw std::runtime_error(
"Mesh format not supported");
599 for (
int c = 0; c < (int)M.facets.nb(); ++c)
601 for (
int lv = 0; lv < F.cols(); ++lv)
603 M.facets.set_vertex(c, lv, F(c, lv));
640 V.resize(M.vertices.nb(), 3);
641 for (
int i = 0; i < (int)M.vertices.nb(); ++i)
643 GEO::vec3 p = M.vertices.point(i);
644 V.row(i) << p[0], p[1], p[2];
646 assert(M.facets.are_simplices());
647 F.resize(M.facets.nb(), 3);
648 for (
int c = 0; c < (int)M.facets.nb(); ++c)
650 for (
int lv = 0; lv < 3; ++lv)
652 F(c, lv) = M.facets.vertex(c, lv);
655 assert(M.cells.are_simplices());
656 T.resize(M.cells.nb(), 4);
657 for (
int c = 0; c < (int)M.cells.nb(); ++c)
659 for (
int lv = 0; lv < 4; ++lv)
661 T(c, lv) = M.cells.vertex(c, lv);
669 const Eigen::MatrixXd &P, Eigen::VectorXd &D)
673 GEO::MeshFacetsAABB aabb_tree(M);
674 compute_unsigned_distance_field(M, aabb_tree, P, D);
675 compute_sign(M, aabb_tree, P, D);
682 assert(F.cols() == 3);
683 assert(
V.cols() == 3);
684 std::array<Eigen::RowVector3d, 4> t;
685 t[3] = Eigen::RowVector3d::Zero(
V.cols());
686 double volume_total = 0;
687 for (
int f = 0; f < F.rows(); ++f)
689 for (
int lv = 0; lv < F.cols(); ++lv)
691 t[lv] =
V.row(F(f, lv));
693 double vol = GEO::Geom::tetra_signed_volume(t[0].data(), t[1].data(), t[2].data(), t[3].data());
696 return -volume_total;
705 for (
int f = 0; f < F.rows(); ++f)
707 F.row(f) = F.row(f).reverse().eval();
717 std::vector<int> vertex_l2g;
718 for (
int c = 0; c < mesh.
n_cells(); ++c)
724 auto poly = std::make_unique<GEO::Mesh>();
727 poly->vertices.create_vertices((triangulated ? nv + nf : nv) + 1);
729 vertex_l2g.reserve(nv);
730 for (
int lf = 0; lf < nf; ++lf)
732 GEO::vector<GEO::index_t> facet_vertices;
736 Eigen::RowVector3d p = mesh.
point(index.vertex);
737 if (vertex_g2l[index.vertex] < 0)
739 vertex_g2l[index.vertex] = vertex_l2g.size();
740 vertex_l2g.push_back(index.vertex);
742 int v1 = vertex_g2l[index.vertex];
743 facet_vertices.push_back(v1);
744 poly->vertices.point(v1) = GEO::vec3(p.data());
749 GEO::vec3 p(0, 0, 0);
750 for (GEO::index_t lv = 0; lv < facet_vertices.size(); ++lv)
752 p += poly->vertices.point(facet_vertices[lv]);
754 p /= facet_vertices.size();
755 int v0 = vertex_l2g.size();
756 vertex_l2g.push_back(0);
757 poly->vertices.point(v0) = p;
758 for (GEO::index_t lv = 0; lv < facet_vertices.size(); ++lv)
760 int v1 = facet_vertices[lv];
761 int v2 = facet_vertices[(lv + 1) % facet_vertices.size()];
762 poly->facets.create_triangle(v0, v1, v2);
767 poly->facets.create_polygon(facet_vertices);
771 Eigen::RowVector3d p = mesh.
kernel(c);
772 poly->vertices.point(nv) = GEO::vec3(p.data());
774 assert(vertex_l2g.size() ==
size_t(triangulated ? nv + nf : nv));
776 for (
int v : vertex_l2g)
781 poly->facets.compute_borders();
782 poly->facets.connect();
784 polys.emplace_back(std::move(poly));
814 M.vertices.create_vertices((
int)mesh.
n_vertices());
815 for (
int i = 0; i < (int)M.vertices.nb(); ++i)
817 auto pt = mesh.
point(i);
818 GEO::vec3 &p = M.vertices.point(i);
824 for (
int f = 0, lf = 0; f < mesh.
n_faces(); ++f)
829 M.facets.create_polygon(nv);
830 for (
int lv = 0; lv < nv; ++lv)
832 M.facets.set_vertex(lf, lv, mesh.
face_vertex(f, lv));
838 typedef std::array<int, 8> Vector8i;
839 Vector8i g2p = {{0, 4, 1, 5, 3, 7, 2, 6}};
840 for (
int c = 0; c < mesh.
n_cells(); ++c)
846 for (
size_t k = 0; k < 8; ++k)
848 lvg[k] = lvp[g2p[k]];
850 std::reverse(lvg.begin(), lvg.end());
852 lvg[0], lvg[1], lvg[2], lvg[3],
853 lvg[4], lvg[5], lvg[6], lvg[7]);
863 GEO::mesh_reorient(M);
869 const Eigen::RowVector3d &kernel, Eigen::MatrixXd &OV, Eigen::MatrixXi &OF, Eigen::MatrixXi &OT)
871 assert(
V.cols() == 3);
872 OV.resize(
V.rows() + 1,
V.cols());
873 OV.topRows(
V.rows()) =
V;
874 OV.bottomRows(1) = kernel;
876 OT.resize(OF.rows(), 4);
877 OT.col(0).setConstant(
V.rows());
884 Eigen::MatrixXd &P, Eigen::MatrixXd *N,
int num_lloyd,
int num_newton)
886 assert(num_samples > 3);
889 GEO::CentroidalVoronoiTesselation CVT(&M);
891 bool was_quiet = GEO::Logger::instance()->is_quiet();
892 GEO::Logger::instance()->set_quiet(
true);
893 CVT.compute_initial_sampling(num_samples);
894 GEO::Logger::instance()->set_quiet(was_quiet);
898 CVT.Lloyd_iterations(num_lloyd);
903 CVT.Newton_iterations(num_newton);
906 P.resize(3, num_samples);
907 std::copy_n(CVT.embedding(0), 3 * num_samples, P.data());
908 P.transposeInPlace();
912 GEO::MeshFacetsAABB aabb(M);
914 for (
int i = 0; i < num_samples; ++i)
916 GEO::vec3 p(P(i, 0), P(i, 1), P(i, 2));
917 GEO::vec3 nearest_point;
919 auto f = aabb.nearest_facet(p, nearest_point, sq_dist);
920 GEO::vec3 n = normalize(GEO::Geom::mesh_facet_normal(M, f));
921 N->row(i) << n[0], n[1], n[2];
931 bool approx_aligned(
const double *a_,
const double *b_,
const double *p_,
const double *q_,
double tol = 1e-6)
934 vec3 a(a_), b(b_), p(p_), q(q_);
935 double da = std::sqrt(Geom::point_segment_squared_distance(a, p, q));
936 double db = std::sqrt(Geom::point_segment_squared_distance(b, p, q));
937 double cos_theta = Geom::cos_angle(b - a, p - q);
938 return (
da < tol && db < tol && std::abs(std::abs(cos_theta) - 1.0) < tol);
946 const Eigen::MatrixXd &BV,
const Eigen::MatrixXi &BE, Eigen::MatrixXi &OE)
948 assert(IV.cols() == 2 || IV.cols() == 3);
949 assert(BV.cols() == 2 || BV.cols() == 3);
950 typedef std::pair<int, int>
Edge;
951 std::vector<Edge> selected;
952 for (
int e1 = 0; e1 < IE.rows(); ++e1)
954 Eigen::RowVector3d a;
956 a.head(IV.cols()) = IV.row(IE(e1, 0));
957 Eigen::RowVector3d b;
959 b.head(IV.cols()) = IV.row(IE(e1, 1));
960 for (
int e2 = 0; e2 < BE.rows(); ++e2)
962 Eigen::RowVector3d p;
964 p.head(BV.cols()) = BV.row(BE(e2, 0));
965 Eigen::RowVector3d q;
967 q.head(BV.cols()) = BV.row(BE(e2, 1));
968 if (approx_aligned(a.data(), b.data(), p.data(), q.data()))
970 selected.emplace_back(IE(e1, 0), IE(e1, 1));
976 OE.resize(selected.size(), 2);
977 for (
int e = 0; e < OE.rows(); ++e)
979 OE.row(e) << selected[e].first, selected[e].second;
986 const Eigen::MatrixXd &vertices,
987 const Eigen::MatrixXi &codim_edges,
988 const Eigen::MatrixXi &
faces,
989 Eigen::VectorXi &codim_vertices)
991 std::vector<bool> is_vertex_codim(vertices.rows(),
true);
992 for (
int i = 0; i < codim_edges.rows(); i++)
994 for (
int j = 0; j < codim_edges.cols(); j++)
996 is_vertex_codim[codim_edges(i, j)] =
false;
999 for (
int i = 0; i <
faces.rows(); i++)
1001 for (
int j = 0; j <
faces.cols(); j++)
1003 is_vertex_codim[
faces(i, j)] =
false;
1006 const auto n_codim_vertices = std::count(is_vertex_codim.begin(), is_vertex_codim.end(),
true);
1007 codim_vertices.resize(n_codim_vertices);
1008 for (
int i = 0, ci = 0; i < vertices.rows(); i++)
1010 if (is_vertex_codim[i])
1012 codim_vertices[ci++] = i;
1018 const Eigen::MatrixXi &tets,
1019 Eigen::MatrixXi &
faces)
1021 std::unordered_set<Eigen::Vector3i, HashMatrix> tri_to_tet(4 * tets.rows());
1022 for (
int i = 0; i < tets.rows(); i++)
1024 tri_to_tet.emplace(tets(i, 0), tets(i, 2), tets(i, 1));
1025 tri_to_tet.emplace(tets(i, 0), tets(i, 3), tets(i, 2));
1026 tri_to_tet.emplace(tets(i, 0), tets(i, 1), tets(i, 3));
1027 tri_to_tet.emplace(tets(i, 1), tets(i, 2), tets(i, 3));
1030 std::vector<Eigen::RowVector3i> faces_vector;
1031 for (
const auto &tri : tri_to_tet)
1034 bool is_surface_triangle =
1035 tri_to_tet.find(Eigen::Vector3i(tri[2], tri[1], tri[0])) == tri_to_tet.end()
1036 && tri_to_tet.find(Eigen::Vector3i(tri[1], tri[0], tri[2])) == tri_to_tet.end()
1037 && tri_to_tet.find(Eigen::Vector3i(tri[0], tri[2], tri[1])) == tri_to_tet.end();
1038 if (is_surface_triangle)
1040 faces_vector.emplace_back(tri[0], tri[1], tri[2]);
1044 faces.resize(faces_vector.size(), 3);
1045 for (
int i = 0; i <
faces.rows(); i++)
1047 faces.row(i) = faces_vector[i];
1052 const Eigen::MatrixXd &vertices,
1053 const Eigen::MatrixXi &tets,
1054 Eigen::MatrixXd &surface_vertices,
1055 Eigen::MatrixXi &tris)
1057 Eigen::MatrixXi full_tris;
1060 std::unordered_map<int, int> full_to_surface;
1061 std::vector<size_t> surface_to_full;
1062 for (
int i = 0; i < full_tris.rows(); i++)
1064 for (
int j = 0; j < full_tris.cols(); j++)
1066 if (full_to_surface.find(full_tris(i, j)) == full_to_surface.end())
1068 full_to_surface[full_tris(i, j)] = surface_to_full.size();
1069 surface_to_full.push_back(full_tris(i, j));
1074 surface_vertices.resize(surface_to_full.size(), 3);
1075 for (
int i = 0; i < surface_to_full.size(); i++)
1077 surface_vertices.row(i) = vertices.row(surface_to_full[i]);
1080 tris.resize(full_tris.rows(), full_tris.cols());
1081 for (
int i = 0; i < tris.rows(); i++)
1083 for (
int j = 0; j < tris.cols(); j++)
1085 tris(i, j) = full_to_surface[full_tris(i, j)];
1091 const std::string &mesh_path,
1092 Eigen::MatrixXd &vertices,
1093 Eigen::VectorXi &codim_vertices,
1094 Eigen::MatrixXi &codim_edges,
1095 Eigen::MatrixXi &
faces)
1097 vertices.resize(0, 0);
1098 codim_vertices.resize(0);
1099 codim_edges.resize(0, 0);
1102 std::string lowername = mesh_path;
1104 lowername.begin(), lowername.end(), lowername.begin(), ::tolower);
1108 Eigen::MatrixXi cells;
1109 std::vector<std::vector<int>> elements;
1110 std::vector<std::vector<double>> weights;
1111 std::vector<int> body_ids;
1112 if (!
MshReader::load(mesh_path, vertices, cells, elements, weights, body_ids))
1114 logger().error(
"Unable to load mesh: {}", mesh_path);
1118 if (cells.cols() == 1)
1119 codim_vertices = cells;
1120 else if (cells.cols() == 2)
1121 codim_edges = cells;
1122 else if (cells.cols() == 3)
1124 else if (cells.cols() == 4)
1126 if (vertices.cols() == 2)
1128 logger().error(
"read_surface_mesh not implemented for 2D quad meshes");
1134 assert(vertices.cols() == 3);
1135 Eigen::MatrixXd surface_vertices;
1137 vertices = surface_vertices;
1142 logger().error(
"read_surface_mesh not implemented for hexahedral and polygonal/polyhedral meshes");
1150 logger().error(
"Unable to load mesh: {}", mesh_path);
1154 else if (!igl::read_triangle_mesh(mesh_path, vertices,
faces))
1157 if (!GEO::mesh_load(mesh_path, mesh))
1159 logger().error(
"Unable to load mesh: {}", mesh_path);
1164 vertices.resize(mesh.vertices.nb(), dim);
1165 for (
int vi = 0; vi < mesh.vertices.nb(); vi++)
1167 const auto &v = mesh.vertices.point(vi);
1168 for (
int vj = 0; vj < dim; vj++)
1170 vertices(vi, vj) = v[vj];
1175 assert(mesh.facets.nb());
1176 int face_cols = mesh.facets.nb_vertices(0);
1177 faces.resize(mesh.facets.nb(), face_cols);
1178 for (
int fi = 0; fi < mesh.facets.nb(); fi++)
1180 assert(face_cols == mesh.facets.nb_vertices(fi));
1181 for (
int fj = 0; fj < mesh.facets.nb_vertices(fi); fj++)
1183 faces(fi, fj) = mesh.facets.vertex(fi, fj);
1195 std::unordered_set<std::vector<int>,
HashVector> boundaries;
1197 auto insert = [&](std::vector<int> v) {
1198 std::sort(v.begin(), v.end());
1199 boundaries.insert(v);
1202 for (
int i = 0; i < cells.rows(); i++)
1204 const auto &cell = cells.row(i);
1205 if (cells.cols() == 3)
1207 insert({{cell(0), cell(1)}});
1208 insert({{cell(1), cell(2)}});
1209 insert({{cell(2), cell(0)}});
1211 else if (cells.cols() == 4 && dim == 2)
1213 insert({{cell(0), cell(1)}});
1214 insert({{cell(1), cell(2)}});
1215 insert({{cell(2), cell(3)}});
1216 insert({{cell(3), cell(0)}});
1218 else if (cells.cols() == 4 && dim == 3)
1220 insert({{cell(0), cell(2), cell(1)}});
1221 insert({{cell(0), cell(3), cell(2)}});
1222 insert({{cell(0), cell(1), cell(3)}});
1223 insert({{cell(1), cell(2), cell(3)}});
1225 else if (cells.cols() == 8)
1227 insert({{cell(0), cell(1), cell(2), cell(3)}});
1228 insert({{cell(1), cell(5), cell(6), cell(3)}});
1229 insert({{cell(5), cell(4), cell(7), cell(6)}});
1230 insert({{cell(0), cell(4), cell(7), cell(3)}});
1231 insert({{cell(0), cell(4), cell(5), cell(1)}});
1232 insert({{cell(2), cell(6), cell(7), cell(3)}});
1236 throw std::runtime_error(
"count_boundary_elements not implemented for polygons");
1240 return boundaries.size();
1245 using namespace GEO;
1246 typedef std::pair<index_t, index_t>
Edge;
1248 if (M.edges.nb() > 0)
1253 if (M.cells.nb() != 0 && M.facets.nb() == 0)
1255 M.cells.compute_borders();
1259 std::vector<std::pair<Edge, index_t>> e2c;
1260 for (index_t f = 0; f < M.facets.nb(); ++f)
1262 for (index_t c = M.facets.corners_begin(f); c < M.facets.corners_end(f); ++c)
1264 index_t v = M.facet_corners.vertex(c);
1265 index_t c2 = M.facets.next_corner_around_facet(f, c);
1266 index_t v2 = M.facet_corners.vertex(c2);
1267 e2c.emplace_back(std::make_pair(std::min(v, v2), std::max(v, v2)), c);
1270 std::sort(e2c.begin(), e2c.end());
1274 Edge prev_e(-1, -1);
1275 for (
const auto &kv : e2c)
1280 M.edges.create_edge(e.first, e.second);
void find_codim_vertices(const Eigen::MatrixXd &vertices, const Eigen::MatrixXi &codim_edges, const Eigen::MatrixXi &faces, Eigen::VectorXi &codim_vertices)
void find_triangle_surface_from_tets(const Eigen::MatrixXi &tets, Eigen::MatrixXi &faces)
std::vector< Eigen::VectorXi > faces
static bool load(const std::string &path, Eigen::MatrixXd &vertices, Eigen::MatrixXi &cells, std::vector< std::vector< int > > &elements, std::vector< std::vector< double > > &weights, std::vector< int > &body_ids)
static bool read(const std::string obj_file_name, std::vector< std::vector< double > > &V, std::vector< std::vector< double > > &TC, std::vector< std::vector< double > > &N, std::vector< std::vector< int > > &F, std::vector< std::vector< int > > &FTC, std::vector< std::vector< int > > &FN, std::vector< std::vector< int > > &L)
Read a mesh from an ascii obj file.
virtual Navigation3D::Index get_index_from_element(int hi, int lf, int lv) const =0
virtual RowVectorNd kernel(const int cell_id) const =0
std::array< int, 8 > get_ordered_vertices_from_hex(const int element_index) const
virtual int n_cell_faces(const int c_id) const =0
virtual Navigation3D::Index next_around_face(Navigation3D::Index idx) const =0
virtual int n_vertices() const =0
number of vertices
bool is_polytope(const int el_id) const
checks if element is polygon compatible
bool is_cube(const int el_id) const
checks if element is cube compatible
virtual RowVectorNd point(const int global_index) const =0
point coordinates
virtual bool is_boundary_face(const int face_global_id) const =0
is face boundary
virtual int n_cells() const =0
number of cells
virtual int n_faces() const =0
number of faces
virtual int n_face_vertices(const int f_id) const =0
number of vertices of a face
virtual int n_cell_vertices(const int c_id) const =0
number of vertices of a cell
virtual int face_vertex(const int f_id, const int lv_id) const =0
id of the face vertex
Eigen::ArrayXd P(const int m, const int p, const Eigen::ArrayXd &z)
bool is_planar(const GEO::Mesh &M, const double tol=1e-5)
Determine if the given mesh is planar (2D or tiny z-range).
void sample_surface(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, int num_samples, Eigen::MatrixXd &P, Eigen::MatrixXd *N=nullptr, int num_lloyd=10, int num_newton=10)
Samples points on a surface.
void orient_closed_surface(const Eigen::MatrixXd &V, Eigen::MatrixXi &F, bool positive=true)
Orient a triangulated surface to have positive volume.
double signed_volume(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F)
Compute the signed volume of a surface mesh.
void reorder_mesh(Eigen::MatrixXd &V, Eigen::MatrixXi &F, const Eigen::VectorXi &C, Eigen::VectorXi &R)
Reorder vertices of a mesh using color tags, so that vertices are ordered by increasing colors.
void orient_normals_2d(GEO::Mesh &M)
Orient facets of a 2D mesh so that each connected component has positive volume.
void extract_parent_edges(const Eigen::MatrixXd &IV, const Eigen::MatrixXi &IE, const Eigen::MatrixXd &BV, const Eigen::MatrixXi &BE, Eigen::MatrixXi &OE)
Extract a set of edges that are overlap with a set given set of parent edges, using vertices position...
ElementType
Type of Element, check [Poly-Spline Finite Element Method] for a complete description.
void extract_triangle_surface_from_tets(const Eigen::MatrixXd &vertices, const Eigen::MatrixXi &tets, Eigen::MatrixXd &surface_vertices, Eigen::MatrixXi &tris)
Extract triangular surface from a tetmesh.
GEO::vec3 mesh_vertex(const GEO::Mesh &M, GEO::index_t v)
Retrieve a 3D vector with the position of a given vertex.
GEO::index_t mesh_create_vertex(GEO::Mesh &M, const GEO::vec3 &p)
void to_geogram_mesh(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, GEO::Mesh &M)
Converts a triangle mesh to a Geogram mesh.
void signed_squared_distances(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const Eigen::MatrixXd &P, Eigen::VectorXd &D)
Computes the signed squared distance from a list of points to a triangle mesh.
GEO::vec3 facet_barycenter(const GEO::Mesh &M, GEO::index_t f)
void from_geogram_mesh(const GEO::Mesh &M, Eigen::MatrixXd &V, Eigen::MatrixXi &F, Eigen::MatrixXi &T)
Extract simplices from a Geogram mesh.
void tertrahedralize_star_shaped_surface(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const Eigen::RowVector3d &kernel, Eigen::MatrixXd &OV, Eigen::MatrixXi &OF, Eigen::MatrixXi &OT)
Tetrahedralize a star-shaped mesh, with a given point in its kernel.
void generate_edges(GEO::Mesh &M)
assing edges to M
bool read_surface_mesh(const std::string &mesh_path, Eigen::MatrixXd &vertices, Eigen::VectorXi &codim_vertices, Eigen::MatrixXi &codim_edges, Eigen::MatrixXi &faces)
read a surface mesh
void extract_polyhedra(const Mesh3D &mesh, std::vector< std::unique_ptr< GEO::Mesh > > &polys, bool triangulated=false)
Extract polyhedra from a 3D volumetric mesh.
int count_faces(const int dim, const Eigen::MatrixXi &cells)
Count the number of boundary elements (triangles for tetmesh and edges for triangle mesh)
void compute_element_tags(const GEO::Mesh &M, std::vector< ElementType > &element_tags)
Compute the type of each facet in a surface mesh.
bool endswith(const std::string &str, const std::string &suffix)
spdlog::logger & logger()
Retrieves the current logger.