diff --git a/include/geometrycentral/surface/dual_mesh.h b/include/geometrycentral/surface/dual_mesh.h new file mode 100644 index 00000000..40c9e7a2 --- /dev/null +++ b/include/geometrycentral/surface/dual_mesh.h @@ -0,0 +1,82 @@ +#pragma once + + +#include "geometrycentral/surface/manifold_surface_mesh.h" +#include "geometrycentral/surface/vertex_position_geometry.h" + + +namespace geometrycentral +{ +namespace surface +{ + +/** + * @brief Builds the dual manifold mesh of the input manifold mesh. + * + * @details This function computes the dual manifold mesh of the given primal + * manifold mesh.\n + * Each face of the primal mesh corresponds to a vertex in the dual mesh. + * If two faces are adjacent (i.e., share an edge) in the primal mesh, + * the corresponding dual vertices are connected by an edge in the dual mesh.\n + * As a consequence, each vertex in the primal mesh corresponds to a face + * in the dual mesh (i.e., the dual vertices corresponding to fan of primal + * faces around the primal vertex).\n + * If the primal manifold mesh has boundaries, it is impossible to preserve + * the vertex-face bijection in both directions.\n + * If the parameter keepBoundaries is set to false, the bijection between + * primal faces and dual vertices is preserved, but the primal vertices at the + * boundary will not have corresponding faces, and the 3D embeddings of the + * primal and dual mesh would have different boundaries.\n + * If the parameter keepBoundaries is set to true, additional dual vertices will + * be created at the primal boundary vertices and at the midpoints of the + * primal boundary edges. This preserves the bijection between primal vertices + * and dual faces, but the dual vertices at boundary will not have corresponding + * primal faces. + * + * @param mesh The mesh whose the dual must be computed. + * @param keepBoundaries Whether or not to create additional dual vertices to preserve the boundary embedding. + * @return std::unique_ptr The resulting dual mesh. + */ +std::unique_ptr +dual_mesh(geometrycentral::surface::ManifoldSurfaceMesh& mesh, + bool keepBoundaries); + +/** + * @brief Builds the dual manifold mesh of the input manifold mesh, with the corresponding + * vertex geometry. + * + * @details This function computes the dual manifold mesh of the given primal + * manifold mesh, with the corresponding vertex geometry.\n + * Each face of the primal mesh corresponds to a vertex in the dual mesh, whose + * position in 3D space is the barycenter of the primal face. + * If two faces are adjacent (i.e., share an edge) in the primal mesh, + * the corresponding dual vertices are connected by an edge in the dual mesh.\n + * As a consequence, each vertex in the primal mesh corresponds to a face + * in the dual mesh (i.e., the dual vertices corresponding to fan of primal + * faces around the primal vertex).\n + * If the primal manifold mesh has boundaries, it is impossible to preserve + * the vertex-face bijection in both directions.\n + * If the parameter keepBoundaries is set to false, the bijection between + * primal faces and dual vertices is preserved, but the primal vertices at the + * boundary will not have corresponding faces, and the 3D embeddings of the + * primal and dual mesh would have different boundaries.\n + * If the parameter keepBoundaries is set to true, additional dual vertices will + * be created at the primal boundary vertices and at the midpoints of the + * primal boundary edges. This preserves the bijection between primal vertices + * and dual faces, but the dual vertices at boundary will not have corresponding + * primal faces. + * + * @param mesh The mesh whose the dual must be computed. + * @param geometry The geometry of the primal mesh. + * @param keepBoundaries Whether or not to create additional dual vertices to preserve the boundary embedding. + * @return std::tuple, + * std::unique_ptr> The tuple containing the dual mesh and its embedded geometry. + */ +std::tuple, + std::unique_ptr> +dual_mesh(geometrycentral::surface::ManifoldSurfaceMesh& mesh, + geometrycentral::surface::VertexPositionGeometry& geometry, + bool keepBoundaries); + +} // namespace surface +} // namespace geometrycentral \ No newline at end of file diff --git a/src/surface/dual_mesh.cpp b/src/surface/dual_mesh.cpp new file mode 100644 index 00000000..0b3eb222 --- /dev/null +++ b/src/surface/dual_mesh.cpp @@ -0,0 +1,279 @@ +#include "geometrycentral/surface/dual_mesh.h" +#include "geometrycentral/surface/surface_mesh_factories.h" + +using namespace geometrycentral; +using namespace geometrycentral::surface; + + +std::tuple, std::unique_ptr> +dual_mesh_nokeep(ManifoldSurfaceMesh& mesh, + VertexPositionGeometry& geometry) +{ + std::vector Barycs; + Barycs.resize(mesh.nFaces()); + for (Face f : mesh.faces()) + { + Barycs[f.getIndex()] = { 0.0, 0.0, 0.0 }; + for (Vertex vf : f.adjacentVertices()) + Barycs[f.getIndex()] += geometry.vertexPositions[vf]; + Barycs[f.getIndex()] /= f.degree(); + } + + std::vector> VFaces; + VFaces.reserve(mesh.nVertices()); + for (Vertex v : mesh.vertices()) + { + std::vector Poly; + Poly.reserve(v.faceDegree() + (v.isBoundary() ? 3 : 0)); + for (Face fv : v.adjacentFaces()) + Poly.emplace_back(fv.getIndex()); + if (v.isBoundary()) + continue; + std::reverse(Poly.begin(), Poly.end()); + VFaces.emplace_back(Poly); + } + + return makeManifoldSurfaceMeshAndGeometry(VFaces, Barycs); +} + +std::tuple, std::unique_ptr> +dual_mesh_keep(ManifoldSurfaceMesh& mesh, + VertexPositionGeometry& geometry) +{ + size_t ExtraVerts = 0; + if (mesh.nBoundaryLoops() > 0) + { + for (BoundaryLoop bl : mesh.boundaryLoops()) + { + for (Vertex v : bl.adjacentVertices()) + ExtraVerts++; + for (Edge e : bl.adjacentEdges()) + ExtraVerts++; + } + } + + std::vector Barycs; + Barycs.reserve(mesh.nFaces() + ExtraVerts); + Barycs.resize(mesh.nFaces()); + for (Face f : mesh.faces()) + { + Barycs[f.getIndex()] = { 0.0, 0.0, 0.0 }; + for (Vertex vf : f.adjacentVertices()) + Barycs[f.getIndex()] += geometry.vertexPositions[vf]; + Barycs[f.getIndex()] /= f.degree(); + } + size_t CurVLen = Barycs.size(); + + std::vector XVertsMap; + XVertsMap.resize(mesh.nVertices(), -1); + std::vector> VFaces; + VFaces.resize(mesh.nVertices()); + for (Vertex v : mesh.vertices()) + { + auto& Poly = VFaces[v.getIndex()]; + Poly.reserve(v.faceDegree() + (v.isBoundary() ? 3 : 0)); + for (Face fv : v.adjacentFaces()) + Poly.emplace_back(fv.getIndex()); + if (v.isBoundary()) + { + int Offset = 0; + for (Face fv : v.adjacentFaces()) + { + if (fv == v.halfedge().face()) + break; + Offset++; + } + std::rotate(Poly.begin(), Poly.begin() + Offset + 1, Poly.end()); + if (XVertsMap[v.getIndex()] == -1) + { + Barycs.emplace_back(0.5 * (geometry.vertexPositions[v] + geometry.vertexPositions[v.halfedge().next().vertex()])); + Poly.push_back(CurVLen); + XVertsMap[v.getIndex()] = CurVLen; + CurVLen++; + } + else + Poly.push_back(XVertsMap[v.getIndex()]); + Barycs.emplace_back(geometry.vertexPositions[v]); + Poly.push_back(CurVLen); + CurVLen++; + int fvid = Poly[0]; + Face fv = mesh.face(fvid); + for (Halfedge he : fv.adjacentHalfedges()) + { + if (he.next().vertex() != v) + continue; + Vertex ov = he.vertex(); + if (XVertsMap[ov.getIndex()] == -1) + { + Barycs.emplace_back(0.5 * (geometry.vertexPositions[v] + geometry.vertexPositions[ov])); + Poly.push_back(CurVLen); + XVertsMap[ov.getIndex()] = CurVLen; + CurVLen++; + } + else + Poly.push_back(XVertsMap[ov.getIndex()]); + } + } + std::reverse(Poly.begin(), Poly.end()); + } + + return makeManifoldSurfaceMeshAndGeometry(VFaces, Barycs); +} + + + + + + + + + + + + + + + + + + + + + + + + + + +std::unique_ptr +dual_mesh_nokeep(ManifoldSurfaceMesh& mesh) +{ + std::vector> VFaces; + VFaces.reserve(mesh.nVertices()); + for (Vertex v : mesh.vertices()) + { + std::vector Poly; + Poly.reserve(v.faceDegree() + (v.isBoundary() ? 3 : 0)); + for (Face fv : v.adjacentFaces()) + Poly.emplace_back(fv.getIndex()); + if (v.isBoundary()) + continue; + std::reverse(Poly.begin(), Poly.end()); + VFaces.emplace_back(Poly); + } + + std::unique_ptr res(new ManifoldSurfaceMesh(VFaces)); + return std::move(res); +} + +std::unique_ptr +dual_mesh_keep(ManifoldSurfaceMesh& mesh) +{ + size_t ExtraVerts = 0; + if (mesh.nBoundaryLoops() > 0) + { + for (BoundaryLoop bl : mesh.boundaryLoops()) + { + for (Vertex v : bl.adjacentVertices()) + ExtraVerts++; + for (Edge e : bl.adjacentEdges()) + ExtraVerts++; + } + } + + size_t CurVLen = mesh.nFaces(); + + std::vector XVertsMap; + XVertsMap.resize(mesh.nVertices(), -1); + std::vector> VFaces; + VFaces.resize(mesh.nVertices()); + for (Vertex v : mesh.vertices()) + { + auto& Poly = VFaces[v.getIndex()]; + Poly.reserve(v.faceDegree() + (v.isBoundary() ? 3 : 0)); + for (Face fv : v.adjacentFaces()) + Poly.emplace_back(fv.getIndex()); + if (v.isBoundary()) + { + int Offset = 0; + for (Face fv : v.adjacentFaces()) + { + if (fv == v.halfedge().face()) + break; + Offset++; + } + std::rotate(Poly.begin(), Poly.begin() + Offset + 1, Poly.end()); + if (XVertsMap[v.getIndex()] == -1) + { + Poly.push_back(CurVLen); + XVertsMap[v.getIndex()] = CurVLen; + CurVLen++; + } + else + Poly.push_back(XVertsMap[v.getIndex()]); + Poly.push_back(CurVLen); + CurVLen++; + int fvid = Poly[0]; + Face fv = mesh.face(fvid); + for (Halfedge he : fv.adjacentHalfedges()) + { + if (he.next().vertex() != v) + continue; + Vertex ov = he.vertex(); + if (XVertsMap[ov.getIndex()] == -1) + { + Poly.push_back(CurVLen); + XVertsMap[ov.getIndex()] = CurVLen; + CurVLen++; + } + else + Poly.push_back(XVertsMap[ov.getIndex()]); + } + } + std::reverse(Poly.begin(), Poly.end()); + } + + std::unique_ptr res(new ManifoldSurfaceMesh(VFaces)); + return std::move(res); +} + + + + + + + + + + + + + + + + + + + + + +std::unique_ptr +geometrycentral::surface::dual_mesh(ManifoldSurfaceMesh& mesh, + bool keepBoundaries) +{ + if (keepBoundaries) + return dual_mesh_keep(mesh); + else + return dual_mesh_nokeep(mesh); +} + +std::tuple, std::unique_ptr> +geometrycentral::surface::dual_mesh(ManifoldSurfaceMesh& mesh, + VertexPositionGeometry& geometry, + bool keepBoundaries) +{ + if (keepBoundaries) + return dual_mesh_keep(mesh, geometry); + else + return dual_mesh_nokeep(mesh, geometry); +} \ No newline at end of file