Dual mesh and dual geometry.

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<geometrycentral::surface::ManifoldSurfaceMesh> The resulting dual mesh.
+ */
+std::unique_ptr<geometrycentral::surface::ManifoldSurfaceMesh>
+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<geometrycentral::surface::ManifoldSurfaceMesh>,
+ * std::unique_ptr<geometrycentral::surface::VertexPositionGeometry>> The tuple containing the dual mesh and its embedded geometry.
+ */
+std::tuple<std::unique_ptr<geometrycentral::surface::ManifoldSurfaceMesh>, 
+           std::unique_ptr<geometrycentral::surface::VertexPositionGeometry>>
+dual_mesh(geometrycentral::surface::ManifoldSurfaceMesh& mesh,
+          geometrycentral::surface::VertexPositionGeometry& geometry,
+          bool keepBoundaries);
+
+} // namespace surface
+} // namespace geometrycentral

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<ManifoldSurfaceMesh>, std::unique_ptr<VertexPositionGeometry>>
+dual_mesh_nokeep(ManifoldSurfaceMesh& mesh,
+                 VertexPositionGeometry& geometry)
+{
+    std::vector<Vector3> 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<std::vector<size_t>> VFaces;
+    VFaces.reserve(mesh.nVertices());
+    for (Vertex v : mesh.vertices())
+    {
+        std::vector<size_t> 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<ManifoldSurfaceMesh>, std::unique_ptr<VertexPositionGeometry>>
+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<Vector3> 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<int> XVertsMap;
+    XVertsMap.resize(mesh.nVertices(), -1);
+    std::vector<std::vector<size_t>> 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<ManifoldSurfaceMesh>
+dual_mesh_nokeep(ManifoldSurfaceMesh& mesh)
+{
+    std::vector<std::vector<size_t>> VFaces;
+    VFaces.reserve(mesh.nVertices());
+    for (Vertex v : mesh.vertices())
+    {
+        std::vector<size_t> 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<ManifoldSurfaceMesh> res(new ManifoldSurfaceMesh(VFaces));
+    return std::move(res);
+}
+
+std::unique_ptr<ManifoldSurfaceMesh>
+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<int> XVertsMap;
+    XVertsMap.resize(mesh.nVertices(), -1);
+    std::vector<std::vector<size_t>> 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<ManifoldSurfaceMesh> res(new ManifoldSurfaceMesh(VFaces));
+    return std::move(res);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+std::unique_ptr<ManifoldSurfaceMesh>
+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<ManifoldSurfaceMesh>, std::unique_ptr<VertexPositionGeometry>>
+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);
+}