Merge pull request #137 from MarkGillespie/remeshing

Remeshing

Author: MarkGillespie

docs/docs/media/spot-triangulation-original.png

@@ -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.

docs/docs/media/spot-triangulation-remeshed.png

@@ -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'

docs/docs/surface/algorithms/remeshing.md

@@ -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<bool> repositionableVertices;
+  void setRepositionableVertices(const VertexData<bool>& repositionableVertex, bool defaultValue = true);
+  void clearRepositionableVertices();
+  bool mayRepositionVertex(Vertex v) const;
+
+  EdgeData<bool> splittableEdges;
+  void setSplittableEdges(const EdgeData<bool>& splittableEdge, bool defaultValue = true);
+  void clearSplittableEdges();
+  bool maySplitEdge(Edge e) const;
+
+  EdgeData<bool> collapsibleEdges;
+  void setCollapsibleEdges(const EdgeData<bool>& collapsibleEdge, bool defaultValue = true);
+  void clearCollapsibleEdges();
+  bool mayCollapseEdge(Edge e) const;
+
+  EdgeData<bool> flippableEdges;
+  void setFlippableEdges(const EdgeData<bool>& flippableEdge, bool defaultValue = true);
+  void clearFlippableEdges();
+  bool mayFlipEdge(Edge e) const;
+
+  FaceData<bool> splittableFaces;
+  void setSplittableFaces(const FaceData<bool>& 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<double>& bSplit);
   Vertex splitFace(Face f, Vector3 newVertexPosition);
   Vertex splitFace(Face f, const std::vector<double>& bSplit, Vector3 newVertexPosition);

docs/mkdocs.yml

@@ -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 <deque>
+
+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

include/geometrycentral/surface/mutation_manager.h

@@ -193,6 +193,7 @@ class SurfaceMesh {
   std::list<std::function<void(const std::vector<size_t>&)>> edgePermuteCallbackList;
   std::list<std::function<void(const std::vector<size_t>&)>> halfedgePermuteCallbackList;
   std::list<std::function<void(const std::vector<size_t>&)>> boundaryLoopPermuteCallbackList;
+  std::list<std::function<void(void)>> compressCallbackList;
 
   // Mesh delete callbacks
   // (this unfortunately seems to be necessary; objects which have registered their callbacks above

include/geometrycentral/surface/remeshing.h

@@ -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