diff --git a/docs/docs/media/spot-triangulation-original.png b/docs/docs/media/spot-triangulation-original.png new file mode 100644 index 00000000..9c544601 Binary files /dev/null and b/docs/docs/media/spot-triangulation-original.png differ diff --git a/docs/docs/media/spot-triangulation-remeshed.png b/docs/docs/media/spot-triangulation-remeshed.png new file mode 100644 index 00000000..dd22b7d0 Binary files /dev/null and b/docs/docs/media/spot-triangulation-remeshed.png differ diff --git a/docs/docs/surface/algorithms/remeshing.md b/docs/docs/surface/algorithms/remeshing.md new file mode 100644 index 00000000..8c1e916c --- /dev/null +++ b/docs/docs/surface/algorithms/remeshing.md @@ -0,0 +1,100 @@ + +| Original | Remeshed | +| - | - | +|![finely_triangulated_spot](/media/spot-triangulation-original.png){: style="max-width: 20em; display: block; margin-left: auto; margin-right: auto;"}|![coarsely_triangulated_spot](/media/spot-triangulation-remeshed.png){: style="max-width: 20em; display: block; margin-left: auto; margin-right: auto;"} + +These routines try to improve mesh quality using repeated rounds of [vertex position smoothing](#tangential-vertex-smoothing), [edge flipping](#extrinsic-delaunay-flipping), and [edge splits/collapses](#edge-length-adjustment). + +??? func "`#!cpp void remesh(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, RemeshOptions options = defaultRemeshOptions);`" + + Remesh a mesh using vertex smoothing along with edge splits, collapses, and flips. Options are passed as a [RemeshOptions](#options) object. + +??? func "`#!cpp void remesh(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, MutationManager& mm, RemeshOptions options = defaultRemeshOptions);`" + + Remesh a mesh using vertex smoothing along with edge splits, collapses, and flips. Options are passed as a [RemeshOptions](#options) object. + All mesh mutations are performed through a `MutationManager` object, which can e.g. keep special vertices fixed or run callback functions when certain mesh mutations occur. + +## Tangential Vertex Smoothing +Vertex smoothing moves each vertex towards the average of its neighborhood. This average can be computed in two ways: in *circumentric* smoothing, every vertex is moved towards the area-weighted average of the circumcenters of its neighboring faces (as in [[Chen & Holst 2011]](https://www.sciencedirect.com/science/article/abs/pii/S0045782510003117)), whereas in *Laplacian* smoothing, every vertex is moved towards the average of the positions of its neighboring vertices. + +In either case we perform *tangential smoothing*, meaning that we only move vertices tangentially along the surface. + +??? func "`#!cpp double smoothByCircumcenter(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, double stepSize = 1, RemeshBoundaryCondition bc = RemeshBoundaryCondition::Tangential);`" + + Move each vertex tangentially towards the area-weighted average of its neighboring face circumcenters. + Returns the average amount that each vertex was moved by. + +??? func "`#!cpp double smoothByCircumcenter(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, MutationManager& mm, double stepSize = 1, RemeshBoundaryCondition bc = RemeshBoundaryCondition::Tangential);`" + + Move each vertex tangentially towards the area-weighted average of its neighboring face circumcenters + , using the provided `MutationManager` to move the vertices + Returns the average amount that each vertex was moved by. + +??? func "`#!cpp double smoothByLaplacian(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, double stepSize = 1, RemeshBoundaryCondition bc = RemeshBoundaryCondition::Tangential);`" + + Move each vertex tangentially towards the average of its neighbors + Returns the average amount that each vertex was moved by. +??? func "`#!cpp double smoothByLaplacian(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, MutationManager& mm, double stepSize = 1, RemeshBoundaryCondition bc = RemeshBoundaryCondition::Tangential);`" + + Move each vertex tangentially towards the average of its neighbors, using the provided `MutationManager` to move the vertices + Returns the average amount that each vertex was moved by. + +### Boundary Conditions +The boundary conditions for smoothing can be set to `RemeshBoundaryCondition::Tangential`, `RemeshBoundaryCondition::Fixed`, or `RemeshBoundaryCondition::Free`. Tangential smoothing allows boundary vertices to move along the tangent direction to the boundary curve, fixed smoothing fixes the boundary vertices, and free smoothing allows boundary vertices to move along any direction in the surface's tangent plane. Leaving boundary vertices free allows the mesh to degenerate, so it is not recommended unless you impose additional constraints. + +## Extrinsic Delaunay Flipping +Delaunay flipping performs edge flips to improve triangle quality. + +!!! warning "No guarantees" + + Unlike the [intrinsic Delaunay flipping routines](/surface/intrinsic_triangulations/basics/#high-level-mutators), extrinsic flipping algorithms are not guaranteed to produce a Delaunay mesh. Nonetheless, these edge flips generally improve triangle quality in practice. + +??? func "`#!cpp size_t fixDelaunay(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom);`" + + Try to make all triangles Delaunay using extrinsic edge flips. + Returns the number of flips performed. + +??? func "`#!cpp size_t fixDelaunay(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, MutationManager& mm);`" + Try to make all triangles Delaunay using extrinsic edge flips, using the provided `MutationManager` to flip edges. + Returns the number of flips performed. + +## Edge Length Adjustment +These routines perform edge splits and collapses to drive mesh edge lengths towards a target value. If curvature adaptation is enabled, this target length is made shorter in high-curvature areas, leading to more resolution there. + +??? func "`#!cpp bool adjustEdgeLengths(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, RemeshOptions options = defaultRemeshOptions);`" + + Apply splits and collapses to adjust edge lengths. + + Reads the following options from `options`, with the following default values and meanings: + + * `#!cpp double options.targetEdgeLength = -1`: the target edge length in flat regions. If `targetEdgeLength` is negative, the target length is set relative to the input mesh's mean length. + * `#!cpp double options.curvatureAdaptation = 0`: how much target edge length should vary due to curvature. Set `curvatureAdaptation` to zero if you want edge lengths to be approximately `targetEdgeLength` everywhere. + * `#!cpp double options.minRelativeLength = 0.05`: the minimum possible edge length in the output mesh. Defined relative to `targetEdgeLength`. + +??? func "`#!cpp bool adjustEdgeLengths(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, MutationManager& mm, RemeshOptions options = defaultRemeshOptions);`" + + Apply splits and collapses to adjust edge lengths. + All splits and collapses are performed using the provided `MutationManager`. + + Reads the following options from `options`, with the following default values and meanings: + + * `#!cpp double options.targetEdgeLength = -1`: the target edge length in flat regions. If `targetEdgeLength` is negative, the target length is set relative to the input mesh's mean length. + * `#!cpp double options.curvatureAdaptation = 0`: how much target edge length should vary due to curvature. Set `curvatureAdaptation` to zero if you want edge lengths to be approximately `targetEdgeLength` everywhere. + * `#!cpp double options.minRelativeLength = 0.05`: the minimum possible edge length in the output mesh. Defined relative to `targetEdgeLength`. + +## Helper Types +### Options +Options are passed in to `remesh` via a `RemeshOptions` object. + +| Field | Default value | Meaning | +|----------------------------------------------------|---------------------------------------|------------------------------------------------------------------------------------------------------------------------------------------------------------| +| `#!cpp double targetEdgeLength` | `-1` | the target edge length in flat regions. If `targetEdgeLength` is negative, the target edge length is set to relative the input mesh's mean edge length | +| `#!cpp size_t maxIterations` | `10` | the maximum number of iterations to run for | +| `#!cpp double curvatureAdaptation` | `0` | how much target length should vary due to curvature. Set curvatureAdaptation to 0 if you want edge lengths to be approximately targetEdgeLength everywhere | +| `#!cpp double minRelativeLength` | `0.05` | the minimum possible edge length allowed in the output mesh. Defined relative to targetEdgeLength | +| `#!cpp RemeshSmoothStyle smoothStyle` | `RemeshSmoothStyle::Circumcentric` | the type of vertex smoothing to use (either `RemeshSmoothStyle::Circumcentric` or `RemeshSmoothStyle::Laplacian`) | +| `#!cpp RemeshBoundaryCondition boundaryCondition` | `RemeshBoundaryCondition::Tangential` | the type of motions allowed for boundary vertices (either `RemeshBoundaryCondition::Fixed`, `RemeshBoundaryCondition::Tangential` or `RemeshBoundaryCondition::Free`) | + +!!! warning "'Fixed' boundary may still move slightly" + + Setting `boundaryCondition` to `RemeshBoundaryCondition::Fixed` only fixes boundary vertices during vertex smoothing. Boundary edges may still be split or collapsed during edge length adjustment, which can also cause the boundary to move slightly. To stop all motion of the boundary, you can pass in a `MutationManager` which marks all boundary edges as not splittable or collapsible. diff --git a/docs/mkdocs.yml b/docs/mkdocs.yml index d99ec6ad..653a2000 100644 --- a/docs/mkdocs.yml +++ b/docs/mkdocs.yml @@ -76,15 +76,16 @@ nav: - 'Barycentric Vector': 'surface/utilities/barycentric_vector.md' - Algorithms: - 'Direction Fields' : 'surface/algorithms/direction_fields.md' + - 'Flip Geodesics' : 'surface/algorithms/flip_geodesics.md' - 'Geodesic Distance' : 'surface/algorithms/geodesic_distance.md' - - 'Tracing Geodesic Paths' : 'surface/algorithms/geodesic_paths.md' - - 'Vector Heat Method' : 'surface/algorithms/vector_heat_method.md' - - 'Surface Centers' : 'surface/algorithms/surface_centers.md' - - 'Surface Sampling' : 'surface/algorithms/surface_sampling.md' - 'Geodesic Centroidal Voronoi Tessellations' : 'surface/algorithms/geodesic_voronoi_tessellations.md' - - 'Robust Geometry' : 'surface/algorithms/robust_geometry.md' - - 'Flip Geodesics' : 'surface/algorithms/flip_geodesics.md' - 'Parameterization' : 'surface/algorithms/parameterization.md' + - 'Remeshing' : 'surface/algorithms/remeshing.md' + - 'Robust Geometry' : 'surface/algorithms/robust_geometry.md' + - 'Surface Centers' : 'surface/algorithms/surface_centers.md' + - 'Surface Sampling' : 'surface/algorithms/surface_sampling.md' + - 'Tracing Geodesic Paths' : 'surface/algorithms/geodesic_paths.md' + - 'Vector Heat Method' : 'surface/algorithms/vector_heat_method.md' - Intrinsic Triangulations: - 'Basics' : 'surface/intrinsic_triangulations/basics.md' - 'Common Subdivision' : 'surface/intrinsic_triangulations/common_subdivision.md' diff --git a/include/geometrycentral/surface/mutation_manager.h b/include/geometrycentral/surface/mutation_manager.h index b912b79b..36827467 100644 --- a/include/geometrycentral/surface/mutation_manager.h +++ b/include/geometrycentral/surface/mutation_manager.h @@ -106,6 +106,45 @@ class MutationManager { // Create a connectivity-only mutation manager, which does not touch any geometry by default MutationManager(ManifoldSurfaceMesh& mesh); + // ====================================================== + // ======== Constraints + // ====================================================== + + // Set up callbacks to update edge constraint data following an edge split + // TODO: work out sane policies for updating after other operations? + void setConstraintCallbacks(); + + // Vertices which are allowed to be repositioned. + // (set to an array which holds true if a vertex may be moved, and false if it should stay fixed) + // Note that this array may be uninitialized, in which case all vertices are allowed to be repositioned + // By default, any newly created vertices are repositionable, unless otherwise specified + // Warning: this prevents vertices from being moved by calls to `repositionVertex()`, and also prevents incident edges + // from being collapsed (since this might move the vertex) + VertexData repositionableVertices; + void setRepositionableVertices(const VertexData& repositionableVertex, bool defaultValue = true); + void clearRepositionableVertices(); + bool mayRepositionVertex(Vertex v) const; + + EdgeData splittableEdges; + void setSplittableEdges(const EdgeData& splittableEdge, bool defaultValue = true); + void clearSplittableEdges(); + bool maySplitEdge(Edge e) const; + + EdgeData collapsibleEdges; + void setCollapsibleEdges(const EdgeData& collapsibleEdge, bool defaultValue = true); + void clearCollapsibleEdges(); + bool mayCollapseEdge(Edge e) const; + + EdgeData flippableEdges; + void setFlippableEdges(const EdgeData& flippableEdge, bool defaultValue = true); + void clearFlippableEdges(); + bool mayFlipEdge(Edge e) const; + + FaceData splittableFaces; + void setSplittableFaces(const FaceData& splittableFace, bool defaultValue = true); + void clearSplittableFaces(); + bool maySplitFace(Face f) const; + // ====================================================== // ======== Low-level mutations @@ -117,7 +156,8 @@ class MutationManager { // many of these, since each adds the burden of a corresponding callback policy to update data. // Move a vertex in 3D space - void repositionVertex(Vertex vert, Vector3 offset); + // Returns true if vertex was repositioned (May return false if, e.g., repositionalVertices[vert] is false) + bool repositionVertex(Vertex vert, Vector3 offset); // Flip an edge. bool flipEdge(Edge e); @@ -131,6 +171,7 @@ class MutationManager { // Returns the new halfedge which points away from the new vertex (so he.vertex() is the new vertex), and is the same // direction as e.halfedge() on the original edge. The halfedge direction of the other part of the new split edge is // also preserved. + // Returns Halfedge() if edge was not split (e.g. if splittableEdges[e] is false) Halfedge splitEdge(Edge e, double tSplit); Halfedge splitEdge(Edge e, Vector3 newVertexPosition); Halfedge splitEdge(Edge e, double tSplit, Vector3 newVertexPosition); @@ -142,6 +183,7 @@ class MutationManager { Vertex collapseEdge(Edge e, double tCollapse, Vector3 newVertexPosition); // Split a face (i.e. insert a vertex into the face), and return the new vertex + // Returns Vertex() if face was not split (e.g. if splittableFaces[f] is false) Vertex splitFace(Face f, const std::vector& bSplit); Vertex splitFace(Face f, Vector3 newVertexPosition); Vertex splitFace(Face f, const std::vector& bSplit, Vector3 newVertexPosition); diff --git a/include/geometrycentral/surface/remeshing.h b/include/geometrycentral/surface/remeshing.h new file mode 100644 index 00000000..65d7e7f7 --- /dev/null +++ b/include/geometrycentral/surface/remeshing.h @@ -0,0 +1,61 @@ +#pragma once + +#include "geometrycentral/surface/barycentric_coordinate_helpers.h" +#include "geometrycentral/surface/manifold_surface_mesh.h" +#include "geometrycentral/surface/mutation_manager.h" +#include "geometrycentral/surface/vertex_position_geometry.h" + +#include + +namespace geometrycentral { +namespace surface { + +enum class RemeshBoundaryCondition { Fixed, Tangential, Free }; +enum class RemeshSmoothStyle { Circumcentric, Laplacian }; + +struct RemeshOptions { + double targetEdgeLength = -1; // the target edge length in flat regions. If `targetEdgeLength` is negative, the target + // edge length is set to relative the input mesh's mean edge length + size_t maxIterations = 10; // the maximum number of iterations to run for + double curvatureAdaptation = 0; // how much target length should vary due to curvature. Set curvatureAdaptation + // to 0 if you want lengths to be approximately targetEdgeLength everywhere + double minRelativeLength = 0.05; // the minimum possible edge length allowed in the output mesh. Defined relative to + // targetEdgeLength + RemeshSmoothStyle smoothStyle = RemeshSmoothStyle::Circumcentric; // smoothing function to use + RemeshBoundaryCondition boundaryCondition = + RemeshBoundaryCondition::Tangential; // allowed movement of boundary vertices +}; +extern const RemeshOptions defaultRemeshOptions; + +// Improve mesh using repeated rounds of edge flipping, vertex position smoothing, and edge splits/collapses +void remesh(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, RemeshOptions options = defaultRemeshOptions); +void remesh(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, MutationManager& mm, + RemeshOptions options = defaultRemeshOptions); + +// Try to make all triangles Delaunay +// Returns the number of flips performed +size_t fixDelaunay(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom); +size_t fixDelaunay(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, MutationManager& mm); + +// Average positions of vertices based on surrounding vertex positions +// Returns the average amount each vertex was moved by +double smoothByLaplacian(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, double stepSize = 1, + RemeshBoundaryCondition bc = RemeshBoundaryCondition::Tangential); +double smoothByLaplacian(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, MutationManager& mm, + double stepSize = 1, RemeshBoundaryCondition bc = RemeshBoundaryCondition::Tangential); + +// Average positions of vertices based on surrounding triangle circumenters as in [Chen & Holst 2011] +// Returns the average amount each vertex was moved by +double smoothByCircumcenter(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, double stepSize = 1, + RemeshBoundaryCondition bc = RemeshBoundaryCondition::Tangential); +double smoothByCircumcenter(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, MutationManager& mm, + double stepSize = 1, RemeshBoundaryCondition bc = RemeshBoundaryCondition::Tangential); + +// applies splits and collapses to adjust edge lengths based on the curvature +bool adjustEdgeLengths(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, + RemeshOptions options = defaultRemeshOptions); +bool adjustEdgeLengths(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, MutationManager& mm, + RemeshOptions options = defaultRemeshOptions); + +} // namespace surface +} // namespace geometrycentral diff --git a/include/geometrycentral/surface/surface_mesh.h b/include/geometrycentral/surface/surface_mesh.h index 03dd23a5..42fe6971 100644 --- a/include/geometrycentral/surface/surface_mesh.h +++ b/include/geometrycentral/surface/surface_mesh.h @@ -193,6 +193,7 @@ class SurfaceMesh { std::list&)>> edgePermuteCallbackList; std::list&)>> halfedgePermuteCallbackList; std::list&)>> boundaryLoopPermuteCallbackList; + std::list> compressCallbackList; // Mesh delete callbacks // (this unfortunately seems to be necessary; objects which have registered their callbacks above diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d71b8c1e..c97da2ed 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -43,6 +43,7 @@ SET(SRCS surface/fast_marching_method.cpp surface/uniformize.cpp surface/parameterize.cpp + surface/remeshing.cpp surface/surgery.cpp surface/simple_idt.cpp surface/exact_polyhedral_geodesics.cpp @@ -119,6 +120,7 @@ SET(HEADERS ${INCLUDE_ROOT}/surface/parameterize.h ${INCLUDE_ROOT}/surface/poisson_disk_sampler.h ${INCLUDE_ROOT}/surface/quadric_error_simplification.h + ${INCLUDE_ROOT}/surface/remeshing.h ${INCLUDE_ROOT}/surface/rich_surface_mesh_data.h ${INCLUDE_ROOT}/surface/rich_surface_mesh_data.ipp ${INCLUDE_ROOT}/surface/polygon_soup_mesh.h diff --git a/src/surface/mutation_manager.cpp b/src/surface/mutation_manager.cpp index 11f4ccd8..31b84835 100644 --- a/src/surface/mutation_manager.cpp +++ b/src/surface/mutation_manager.cpp @@ -10,15 +10,93 @@ namespace surface { // ====================================================== MutationManager::MutationManager(ManifoldSurfaceMesh& mesh_, VertexPositionGeometry& geometry_) - : mesh(mesh_), geometry(&geometry_) {} + : mesh(mesh_), geometry(&geometry_) { + setConstraintCallbacks(); +} -MutationManager::MutationManager(ManifoldSurfaceMesh& mesh_) : mesh(mesh_) {} +MutationManager::MutationManager(ManifoldSurfaceMesh& mesh_) : mesh(mesh_) { setConstraintCallbacks(); } + +// ====================================================== +// ======== Constraints +// ====================================================== + +void MutationManager::setConstraintCallbacks() { + // == Register a callback to maintain edge constraint data after edge splits + auto cacheConstraintsBeforeEdgeSplit = [&](Edge oldE, double tSplit) -> std::array { + bool wasSplittable = maySplitEdge(oldE); + bool wasCollapsible = mayCollapseEdge(oldE); + bool wasFlippable = mayFlipEdge(oldE); + return {wasSplittable, wasCollapsible, wasFlippable}; + }; + auto updateConstraintsAfterEdgeSplit = [&](Halfedge newHe1, Halfedge newHe2, double tSplit, + std::array oldConstraints) { + if (splittableEdges.size() > 0) { + splittableEdges[newHe1.edge()] = oldConstraints[0]; + splittableEdges[newHe2.edge()] = oldConstraints[0]; + } + if (collapsibleEdges.size() > 0) { + collapsibleEdges[newHe1.edge()] = oldConstraints[1]; + collapsibleEdges[newHe2.edge()] = oldConstraints[1]; + } + if (flippableEdges.size() > 0) { + flippableEdges[newHe1.edge()] = oldConstraints[2]; + flippableEdges[newHe2.edge()] = oldConstraints[2]; + } + }; + registerEdgeSplitHandlers(cacheConstraintsBeforeEdgeSplit, updateConstraintsAfterEdgeSplit); +} + +void MutationManager::setRepositionableVertices(const VertexData& repositionableVertices_, bool defaultValue) { + repositionableVertices = repositionableVertices_; + repositionableVertices.setDefault(defaultValue); +} +void MutationManager::clearRepositionableVertices() { repositionableVertices = VertexData(); } +bool MutationManager::mayRepositionVertex(Vertex v) const { + return repositionableVertices.size() == 0 || repositionableVertices[v]; +} + +void MutationManager::setSplittableEdges(const EdgeData& splittableEdges_, bool defaultValue) { + splittableEdges = splittableEdges_; + splittableEdges.setDefault(defaultValue); +} +void MutationManager::clearSplittableEdges() { splittableEdges = EdgeData(); } +bool MutationManager::maySplitEdge(Edge e) const { return splittableEdges.size() == 0 || splittableEdges[e]; } + +void MutationManager::setCollapsibleEdges(const EdgeData& collapsibleEdges_, bool defaultValue) { + collapsibleEdges = collapsibleEdges_; + collapsibleEdges.setDefault(defaultValue); +} +void MutationManager::clearCollapsibleEdges() { collapsibleEdges = EdgeData(); } +bool MutationManager::mayCollapseEdge(Edge e) const { + if (collapsibleEdges.size() > 0 && !collapsibleEdges[e]) return false; + if (repositionableVertices.size() > 0) { + for (Vertex v : e.adjacentVertices()) { + if (!repositionableVertices[v]) return false; + } + } + return collapsibleEdges.size() == 0 || collapsibleEdges[e]; +} + +void MutationManager::setFlippableEdges(const EdgeData& flippableEdges_, bool defaultValue) { + flippableEdges = flippableEdges_; + flippableEdges.setDefault(defaultValue); +} +void MutationManager::clearFlippableEdges() { flippableEdges = EdgeData(); } +bool MutationManager::mayFlipEdge(Edge e) const { return flippableEdges.size() == 0 || flippableEdges[e]; } + +void MutationManager::setSplittableFaces(const FaceData& splittableFaces_, bool defaultValue) { + splittableFaces = splittableFaces_; + splittableFaces.setDefault(defaultValue); +} +void MutationManager::clearSplittableFaces() { splittableFaces = FaceData(); } +bool MutationManager::maySplitFace(Face f) const { return splittableFaces.size() == 0 || splittableFaces[f]; } // ====================================================== // ======== Low-level mutations // ====================================================== -void MutationManager::repositionVertex(Vertex vert, Vector3 offset) { +bool MutationManager::repositionVertex(Vertex vert, Vector3 offset) { + if (!mayRepositionVertex(vert)) return false; // Invoke callbacks for (VertexRepositionPolicy* policy : vertexRepositionPolicies) { @@ -26,10 +104,12 @@ void MutationManager::repositionVertex(Vertex vert, Vector3 offset) { } geometry->vertexPositions[vert] += offset; + return true; } // Flip an edge. bool MutationManager::flipEdge(Edge e) { + if (!mayFlipEdge(e)) return false; // First do a test flip to see if its possible // TODO implement canFlip() to avoid this @@ -60,6 +140,8 @@ bool MutationManager::flipEdge(Edge e) { } Halfedge MutationManager::splitEdge(Edge e, double tSplit) { + if (!maySplitEdge(e)) return Halfedge(); + Vector3 newPos{0., 0., 0.}; if (geometry) { VertexData& pos = geometry->vertexPositions; @@ -69,6 +151,7 @@ Halfedge MutationManager::splitEdge(Edge e, double tSplit) { } Halfedge MutationManager::splitEdge(Edge e, Vector3 newVertexPosition) { + if (!maySplitEdge(e)) return Halfedge(); double tSplit = -1; GC_SAFETY_ASSERT(geometry, "must have geometry to split by position"); @@ -84,6 +167,7 @@ Halfedge MutationManager::splitEdge(Edge e, Vector3 newVertexPosition) { } Halfedge MutationManager::splitEdge(Edge e, double tSplit, Vector3 newVertexPosition) { + if (!maySplitEdge(e)) return Halfedge(); // Invoke before callbacks for (EdgeSplitPolicy* policy : edgeSplitPolicies) { @@ -109,6 +193,8 @@ Halfedge MutationManager::splitEdge(Edge e, double tSplit, Vector3 newVertexPosi // Collapse an edge. // Returns the new vertex if the edge could be collapsed, and Vertex() otherwise Vertex MutationManager::collapseEdge(Edge e, double tCollapse) { + if (!mayCollapseEdge(e)) return Vertex(); + Vector3 newPos{0., 0., 0.}; if (geometry) { // Find the nearest tCoord @@ -119,6 +205,7 @@ Vertex MutationManager::collapseEdge(Edge e, double tCollapse) { } Vertex MutationManager::collapseEdge(Edge e, Vector3 newVertexPosition) { + if (!mayCollapseEdge(e)) return Vertex(); double tCollapse = -1; GC_SAFETY_ASSERT(geometry, "must have geometry to split by position"); @@ -133,6 +220,7 @@ Vertex MutationManager::collapseEdge(Edge e, Vector3 newVertexPosition) { return collapseEdge(e, tCollapse, newVertexPosition); } Vertex MutationManager::collapseEdge(Edge e, double tCollapse, Vector3 newVertexPosition) { + if (!mayCollapseEdge(e)) return Vertex(); // Invoke before callbacks // TODO need to handle possiblity that collapse fails -- check before calling @@ -158,6 +246,8 @@ Vertex MutationManager::collapseEdge(Edge e, double tCollapse, Vector3 newVertex // Split a face (i.e. insert a vertex into the face) Vertex MutationManager::splitFace(Face f, const std::vector& bSplit) { + if (!maySplitFace(f)) return Vertex(); + Vector3 newPos = Vector3::zero(); if (geometry) { size_t iV = 0; @@ -172,6 +262,8 @@ Vertex MutationManager::splitFace(Face f, const std::vector& bSplit) { } Vertex MutationManager::splitFace(Face f, Vector3 newVertexPosition) { + if (!maySplitFace(f)) return Vertex(); + // TODO throw std::runtime_error("Face split based on vertex position not implemented yet"); return Vertex(); @@ -179,6 +271,8 @@ Vertex MutationManager::splitFace(Face f, Vector3 newVertexPosition) { Vertex MutationManager::splitFace(Face f, const std::vector& bSplit, Vector3 newVertexPosition) { + if (!maySplitFace(f)) return Vertex(); + // Invoke before callbacks for (FaceSplitPolicy* policy : faceSplitPolicies) { policy->beforeFaceSplit(f, bSplit); diff --git a/src/surface/remeshing.cpp b/src/surface/remeshing.cpp new file mode 100644 index 00000000..fa73e070 --- /dev/null +++ b/src/surface/remeshing.cpp @@ -0,0 +1,429 @@ +#include "geometrycentral/surface/remeshing.h" + +namespace geometrycentral { +namespace surface { + +// The default trace options +const RemeshOptions defaultRemeshOptions; + +void remesh(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, RemeshOptions options) { + MutationManager mm(mesh, geom); + return remesh(mesh, geom, mm, options); +} + +void remesh(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, MutationManager& mm, RemeshOptions options) { + if (options.targetEdgeLength < 0) { + double meanLength = 0; + geom.requireEdgeLengths(); + for (Edge e : mesh.edges()) { + meanLength += geom.edgeLengths[e]; + } + geom.unrequireEdgeLengths(); + meanLength /= mesh.nEdges(); + + options.targetEdgeLength = abs(options.targetEdgeLength * meanLength); + } + + bool doConnectivityChanges = true; + + for (size_t iIt = 0; iIt < options.maxIterations; iIt++) { + if (doConnectivityChanges) { + doConnectivityChanges = adjustEdgeLengths(mesh, geom, mm, options); + } + + size_t nFlips = fixDelaunay(mesh, geom, mm); + double flowDist; + switch (options.smoothStyle) { + case RemeshSmoothStyle::Circumcentric: + flowDist = smoothByCircumcenter(mesh, geom, mm, 1, options.boundaryCondition); + break; + case RemeshSmoothStyle::Laplacian: + flowDist = smoothByLaplacian(mesh, geom, mm, 1, options.boundaryCondition); + break; + } + + // std::cout << iIt << " : " << changedConnectivity << " " << nFlips << " " << flowDist << std::endl; + if ((nFlips == 0) && (flowDist < 0.01)) break; + } + geom.refreshQuantities(); +} + +Vector3 vertexNormal(VertexPositionGeometry& geom, Vertex v, MutationManager& mm) { + Vector3 totalNormal = Vector3::zero(); + for (Corner c : v.adjacentCorners()) { + Vector3 cornerNormal = geom.cornerAngle(c) * geom.faceNormal(c.face()); + totalNormal += cornerNormal; + } + return normalize(totalNormal); +} + +Vector3 boundaryVertexTangent(VertexPositionGeometry& geom, Vertex v, MutationManager& mm) { + if (v.isBoundary()) { + auto edgeVec = [&](Edge e) -> Vector3 { + return (geom.vertexPositions[e.halfedge().tipVertex()] - geom.vertexPositions[e.halfedge().tailVertex()]) + .normalize(); + }; + + Vector3 totalTangent = Vector3::zero(); + for (Edge e : v.adjacentEdges()) { + if (e.isBoundary()) { + totalTangent += edgeVec(e); + } + } + return totalTangent.normalize(); + } else { + return Vector3::zero(); + } +} + +inline Vector3 projectToPlane(Vector3 v, Vector3 norm) { return v - norm * dot(norm, v); } +inline Vector3 projectToLine(Vector3 v, Vector3 tangent) { return tangent * dot(tangent, v); } + +inline Vector3 edgeMidpoint(SurfaceMesh& mesh, VertexPositionGeometry& geom, Edge e) { + Vector3 endPos1 = geom.vertexPositions[e.halfedge().tailVertex()]; + Vector3 endPos2 = geom.vertexPositions[e.halfedge().tipVertex()]; + return (endPos1 + endPos2) / 2; +} + +Vector3 findCircumcenter(Vector3 p1, Vector3 p2, Vector3 p3) { + // barycentric coordinates of circumcenter + double a = (p3 - p2).norm(); + double b = (p3 - p1).norm(); + double c = (p2 - p1).norm(); + double a2 = a * a; + double b2 = b * b; + double c2 = c * c; + Vector3 circumcenterLoc{a2 * (b2 + c2 - a2), b2 * (c2 + a2 - b2), c2 * (a2 + b2 - c2)}; + // normalize to sum of 1 + circumcenterLoc = normalizeBarycentric(circumcenterLoc); + + // change back to space + return circumcenterLoc[0] * p1 + circumcenterLoc[1] * p2 + circumcenterLoc[2] * p3; +} + +Vector3 findCircumcenter(VertexPositionGeometry& geom, Face f) { + // retrieve the face's vertices + int index = 0; + Vector3 p[3]; + for (Vertex v0 : f.adjacentVertices()) { + p[index] = geom.vertexPositions[v0]; + index++; + } + return findCircumcenter(p[0], p[1], p[2]); +} + +// Returns the barycenter for faces incident on a nonflippable edge (e.g. a boundary edge), and the circumcenter for all +// other faces +Vector3 findODTCenter(VertexPositionGeometry& geom, Face f, MutationManager& mm) { + Vector3 p0 = geom.vertexPositions[f.halfedge().tailVertex()]; + Vector3 p1 = geom.vertexPositions[f.halfedge().tipVertex()]; + Vector3 p2 = geom.vertexPositions[f.halfedge().next().tipVertex()]; + + for (Edge e : f.adjacentEdges()) { + if (e.isBoundary() || !mm.mayFlipEdge(e)) { + // e is not flippable. return barycenter + return (p0 + p1 + p2) / 3.; + } + } + return findCircumcenter(p0, p1, p2); +} + +bool isDelaunay(VertexPositionGeometry& geom, Edge e) { + float angle1 = geom.cornerAngle(e.halfedge().next().next().corner()); + float angle2 = geom.cornerAngle(e.halfedge().twin().next().next().corner()); + return angle1 + angle2 <= PI; +} + +inline double diamondAngle(Vector3 a, Vector3 b, Vector3 c, Vector3 d) // dihedral angle at edge a-b +{ + Vector3 n1 = cross(b - a, c - a); + Vector3 n2 = cross(b - d, a - d); + return PI - angle(n1, n2); +} + +inline bool checkFoldover(Vector3 a, Vector3 b, Vector3 c, Vector3 x, double angle) { + return diamondAngle(a, b, c, x) < angle; +} + +bool shouldCollapse(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, Edge e) { + std::vector edgesToCheck; + Vertex v1 = e.halfedge().vertex(); + Vertex v2 = e.halfedge().twin().vertex(); + + // find (halfedge) link around the edge, starting with those surrounding v1 + for (Halfedge he : v1.outgoingHalfedges()) { + if (he.next().tailVertex() != v2 && he.next().tipVertex() != v2) { + edgesToCheck.push_back(he.next()); + } + } + + // link around v2 + for (Halfedge he : v2.outgoingHalfedges()) { + if (he.next().tailVertex() != v1 && he.next().tipVertex() != v1) { + edgesToCheck.push_back(he.next()); + } + } + + // see if the point that would form after a collapse would cause a major foldover with surrounding edges + Vector3 midpoint = edgeMidpoint(mesh, geom, e); + for (Halfedge he0 : edgesToCheck) { + Halfedge heT = he0.twin(); + Vertex v1 = heT.tailVertex(); + Vertex v2 = heT.tipVertex(); + Vertex v3 = heT.next().tipVertex(); + Vector3 a = geom.vertexPositions[v1]; + Vector3 b = geom.vertexPositions[v2]; + Vector3 c = geom.vertexPositions[v3]; + if (checkFoldover(a, b, c, midpoint, 2)) { + return false; + } + } + return true; +} + +// Warning: requires that geom.vertexDualAreas and geom.vertexMeanCurvatures are filled in with accurate data +double getSmoothMeanCurvature(VertexPositionGeometry& geom, Vertex v) { + double A = geom.vertexDualAreas[v]; + double S = geom.vertexMeanCurvatures[v]; + double K = S / A; + return K; +} + +// flatLength: specifies how long the target edge length should be in flat regions +// curvatureAdaptation: controls how much variation in target length occurs due to curvature +double findMeanTargetL(SurfaceMesh& mesh, VertexPositionGeometry& geom, Edge e, double flatLength, + double curvatureAdaptation) { + double averageH = 0; + for (Vertex v : e.adjacentVertices()) { + averageH += getSmoothMeanCurvature(geom, v); + } + averageH /= 2; + double L = flatLength / (fabs(averageH) * curvatureAdaptation + 1); + return L; +} + + +size_t fixDelaunay(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom) { + MutationManager mm(mesh, geom); + return fixDelaunay(mesh, geom, mm); +} + +size_t fixDelaunay(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, MutationManager& mm) { + // Logic duplicated from surface/intrinsic_triangulation.cpp + + std::deque edgesToCheck; // queue of edges to check if Delaunay + EdgeData inQueue(mesh, true); // true if edge is currently in edgesToCheck + + // start with all edges + for (Edge e : mesh.edges()) { + edgesToCheck.push_back(e); + } + + // counter and limit for number of flips + size_t flipMax = 100 * mesh.nVertices(); + size_t nFlips = 0; + while (!edgesToCheck.empty() && nFlips < flipMax) { + Edge e = edgesToCheck.front(); + edgesToCheck.pop_front(); + inQueue[e] = false; + + if (e.isBoundary() || isDelaunay(geom, e)) continue; + + // if not Delaunay, try to flip edge + bool wasFlipped = mm.flipEdge(e); + + if (!wasFlipped) continue; + + nFlips++; + + // Add neighbors to queue, as they may need flipping now + Halfedge he = e.halfedge(); + std::array neighboringEdges{he.next().edge(), he.next().next().edge(), he.twin().next().edge(), + he.twin().next().next().edge()}; + for (Edge nE : neighboringEdges) { + if (!inQueue[nE]) { + edgesToCheck.push_back(nE); + inQueue[nE] = true; + } + } + } + return nFlips; +} + +double smoothByLaplacian(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, double stepSize, + RemeshBoundaryCondition bc) { + MutationManager mm(mesh, geom); + return smoothByLaplacian(mesh, geom, mm, stepSize, bc); +} + +double smoothByLaplacian(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, MutationManager& mm, double stepSize, + RemeshBoundaryCondition bc) { + VertexData vertexOffsets(mesh); + for (Vertex v : mesh.vertices()) { + // calculate average of surrounding vertices + Vector3 avgNeighbor = Vector3::zero(); + for (Vertex j : v.adjacentVertices()) { + avgNeighbor += geom.vertexPositions[j]; + } + avgNeighbor /= v.degree(); + + Vector3 updateDirection = avgNeighbor - geom.vertexPositions[v]; + + // project updateDirection onto space of allowed movements + Vector3 stepDir; + if (v.isBoundary()) { + switch (bc) { + case RemeshBoundaryCondition::Fixed: + stepDir = Vector3::zero(); + break; + case RemeshBoundaryCondition::Tangential: + // for free boundary vertices, project the average to the boundary tangent line + stepDir = projectToLine(updateDirection, boundaryVertexTangent(geom, v, mm)); + break; + case RemeshBoundaryCondition::Free: + // for free boundary vertices, project the average to the surface tangent plane + stepDir = projectToPlane(updateDirection, vertexNormal(geom, v, mm)); + break; + } + } else { + // for interior vertices, project the average to the tangent plane + stepDir = projectToPlane(updateDirection, vertexNormal(geom, v, mm)); + } + vertexOffsets[v] = stepSize * stepDir; + } + + // update final vertices + double totalMovement = 0; + for (Vertex v : mesh.vertices()) { + bool didMove = mm.repositionVertex(v, vertexOffsets[v]); + if (didMove) { + totalMovement += vertexOffsets[v].norm(); + } + } + return totalMovement / mesh.nVertices(); +} + +double smoothByCircumcenter(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, double stepSize, + RemeshBoundaryCondition bc) { + MutationManager mm(mesh, geom); + return smoothByCircumcenter(mesh, geom, mm, stepSize, bc); +} + +double smoothByCircumcenter(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, MutationManager& mm, + double stepSize, RemeshBoundaryCondition bc) { + geom.requireFaceAreas(); + VertexData vertexOffsets(mesh); + for (Vertex v : mesh.vertices()) { + Vector3 updateDirection = Vector3::zero(); + for (Face f : v.adjacentFaces()) { + // add the circumcenter weighted by face area to the update direction + Vector3 circum = findODTCenter(geom, f, mm); + updateDirection += geom.faceArea(f) * (circum - geom.vertexPositions[v]); + } + updateDirection /= (6 * geom.vertexDualArea(v)); + + // project updateDirection onto space of allowed movements + Vector3 stepDir; + if (v.isBoundary()) { + switch (bc) { + case RemeshBoundaryCondition::Fixed: + stepDir = Vector3::zero(); + break; + case RemeshBoundaryCondition::Tangential: + // for free boundary vertices, project the average to the boundary tangent line + stepDir = projectToLine(updateDirection, boundaryVertexTangent(geom, v, mm)); + break; + case RemeshBoundaryCondition::Free: + // for free boundary vertices, project the average to the surface tangent plane + stepDir = projectToPlane(updateDirection, vertexNormal(geom, v, mm)); + break; + } + } else { + // for interior vertices, project the average to the tangent plane + stepDir = projectToPlane(updateDirection, vertexNormal(geom, v, mm)); + } + vertexOffsets[v] = stepSize * stepDir; + } + + // update final vertices + double totalMovement = 0; + for (Vertex v : mesh.vertices()) { + bool didMove = mm.repositionVertex(v, vertexOffsets[v]); + if (didMove) { + totalMovement += vertexOffsets[v].norm(); + } + } + return totalMovement / mesh.nVertices(); +} // namespace surface + + +bool adjustEdgeLengths(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, RemeshOptions options) { + MutationManager mm(mesh, geom); + return adjustEdgeLengths(mesh, geom, mm, options); +} + +bool adjustEdgeLengths(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geom, MutationManager& mm, + RemeshOptions options) { + geom.requireVertexDualAreas(); + geom.requireVertexMeanCurvatures(); + + double curvatureAdaptation = options.curvatureAdaptation; + bool useCurvatureAdaptation = curvatureAdaptation > 0; + double flatLength = options.targetEdgeLength; + double minLength = options.minRelativeLength * options.targetEdgeLength; + + bool didSplitOrCollapse = false; + // queues of edges to CHECK to change + std::vector toSplit; + std::vector toCollapse; + + for (Edge e : mesh.edges()) { + toSplit.push_back(e); + } + + // actually splitting + while (!toSplit.empty()) { + Edge e = toSplit.back(); + toSplit.pop_back(); + double length_e = geom.edgeLength(e); + double threshold = + (useCurvatureAdaptation) ? findMeanTargetL(mesh, geom, e, flatLength, curvatureAdaptation) : flatLength; + if (length_e > minLength && length_e > threshold * 1.5) { + Vector3 newPos = edgeMidpoint(mesh, geom, e); + Halfedge he = mm.splitEdge(e, newPos); + if (he != Halfedge()) { + didSplitOrCollapse = true; + } + } else { + toCollapse.push_back(e); + } + } + + // actually collapsing + while (!toCollapse.empty()) { + Edge e = toCollapse.back(); + toCollapse.pop_back(); + if (e == Edge() || e.isDead()) continue; // make sure it exists + + double threshold = + (useCurvatureAdaptation) ? findMeanTargetL(mesh, geom, e, flatLength, curvatureAdaptation) : flatLength; + if (geom.edgeLength(e) < threshold * 0.5) { + Vector3 newPos = edgeMidpoint(mesh, geom, e); + if (shouldCollapse(mesh, geom, e)) { + Vertex v = mm.collapseEdge(e, newPos); + if (v != Vertex()) { + didSplitOrCollapse = true; + } + } + } + } + + geom.unrequireVertexDualAreas(); + geom.unrequireVertexMeanCurvatures(); + + mesh.compress(); + return didSplitOrCollapse; +} + +} // namespace surface +} // namespace geometrycentral diff --git a/src/surface/surface_mesh.cpp b/src/surface/surface_mesh.cpp index 5f66a220..4d0bdebe 100644 --- a/src/surface/surface_mesh.cpp +++ b/src/surface/surface_mesh.cpp @@ -1940,6 +1940,10 @@ void SurfaceMesh::compress() { compressFaces(); compressVertices(); isCompressedFlag = true; + + for (auto& f : compressCallbackList) { + f(); + } } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index f8ef6d7e..598d6667 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -81,6 +81,12 @@ FetchContent_MakeAvailable(googletest) set_property(TARGET gtest PROPERTY CXX_STANDARD 11) +# silence warning about deprecated copies in gtest +# https://github.com/abseil/abseil-cpp/issues/948#issuecomment-846592426 +if(NOT "${CMAKE_CXX_COMPILER_ID}" STREQUAL "MSVC") + target_compile_options(gtest INTERFACE -Wno-deprecated-copy) +endif() + ### Set (hardcode) a variable for where to find the assets set(REL_ASSETS_DIR "${CMAKE_CURRENT_SOURCE_DIR}/assets/")