Merge pull request #106 from MarkGillespie/IntrinsicEdgeSplits

Improvements to integer coordinate edge splits

include/geometrycentral/surface/integer_coordinates_intrinsic_triangulation.h

@@ -48,6 +48,8 @@ class IntegerCoordinatesIntrinsicTriangulation : public IntrinsicTriangulation {
   EdgeData<std::vector<SurfacePoint>> traceAllInputEdgesAlongIntrinsic() override;
   std::vector<SurfacePoint> traceInputHalfedgeAlongIntrinsic(Halfedge inputHe) override;
 
+  bool checkEdgeOriginal(Edge e) const override;
+
   // ======================================================
   // ======== Low-Level Mutators
   // ======================================================
@@ -64,16 +66,17 @@ class IntegerCoordinatesIntrinsicTriangulation : public IntrinsicTriangulation {
 
   Halfedge splitEdge(Halfedge he, double tSplit) override;
 
-  // Check if an edge can be flipped geometrically, as defined by the (relative) signed areas of the resulting triangles; positive values mean flippable.
+  // Check if an edge can be flipped geometrically, as defined by the (relative) signed areas of the resulting
+  // triangles; positive values mean flippable.
   double checkFlip(Edge e);
 
   // Insert circumcenter or split segment
   Vertex insertCircumcenterOrSplitSegment(Face f, bool verbose = false);
 
   Vertex splitFace(Face f, Vector3 bary, bool verbose = false);
   Vertex splitEdge(Edge e, double bary, bool verbose = false);
-  Vertex splitInteriorEdge(Edge e, double bary, bool verbose = false);
-  Vertex splitBoundaryEdge(Edge e, double bary, bool verbose = false);
+  Halfedge splitInteriorEdge(Halfedge he, double bary, bool verbose = false);
+  Halfedge splitBoundaryEdge(Halfedge he, double bary, bool verbose = false);
 
   // Move a vertex `v` in direction `vec`, represented as a vector in the
   // vertex's tangent space.
@@ -100,6 +103,10 @@ class IntegerCoordinatesIntrinsicTriangulation : public IntrinsicTriangulation {
   // the edges connecting this new point to f's vertices.
   std::pair<SurfacePoint, std::array<int, 3>> computeFaceSplitData(Face f, Vector3 bary, bool verbose = false);
 
+  // Computes the corresponding point on the input mesh, as well as which segment the point belongs to.
+  // The segment is specified by an integer 0 <= segmentIndex <= normalCoordinates[he.edge()]
+  std::pair<SurfacePoint, size_t> computeEdgeSplitData(Halfedge he, double tBary);
+
   // Compute the number of vertices in the common subdivision
   // i.e. intrinsicMesh->nVertices() plus the sum of the normal coordinates
   size_t nSubdividedVertices() const;
@@ -135,10 +142,6 @@ class IntegerCoordinatesIntrinsicTriangulation : public IntrinsicTriangulation {
   void constructCommonSubdivision() override;
 };
 
-// Compute the cotan weight of halfedge ij in terms of the lengths of its
-// neighbors
-double halfedgeCotanWeight(double lij, double ljk, double lki);
-
 FaceData<Vector2> interpolateTangentVectorsB(const IntegerCoordinatesIntrinsicTriangulation& tri,
                                              const CommonSubdivision& cs, const FaceData<Vector2>& dataB);
 

include/geometrycentral/surface/intrinsic_triangulation.h

@@ -66,6 +66,7 @@ class IntrinsicTriangulation : public EdgeLengthGeometry {
   // this possiblity.
   EdgeData<bool> markedEdges;
   void setMarkedEdges(const EdgeData<bool>& markedEdges);
+  void clearMarkedEdges();
   // Is this a marked or boundary edge?
   bool isFixed(Edge e) const;
   bool isOnFixedEdge(Vertex v) const; // boundary vertex or on fixed edge
@@ -128,6 +129,12 @@ class IntrinsicTriangulation : public EdgeLengthGeometry {
   // face Otherwise return Face()
   Face getParentFace(Face f) const;
 
+  // Check if edge is shared with input mesh
+  virtual bool checkEdgeOriginal(Edge e) const = 0;
+
+  // Immediate computation of corner angle
+  double getCornerAngle(Corner c) const;
+
   // ======================================================
   // ======== High-Level Mutators
   // ======================================================
@@ -185,6 +192,13 @@ class IntrinsicTriangulation : public EdgeLengthGeometry {
   // Split an edge
   virtual Halfedge splitEdge(Halfedge he, double tSplit) = 0;
 
+
+  // ==== Misc
+  // Recover t-values after tracing
+  // Note that really we ought to just report these back from the tracing routine itself, which computes them
+  // internally. We don't have a nice API for passing that data around, so this lazily recovers it
+  std::vector<double> recoverTraceTValues(const std::vector<SurfacePoint>& edgeTrace);
+
   // ======================================================
   // ======== Callbacks
   // ======================================================

include/geometrycentral/surface/normal_coordinates.h

@@ -49,7 +49,7 @@ class NormalCoordinates {
 
   // Count of edges emanating from each corner (nonnegative)
   VertexData<int> roundaboutDegrees;
-  
+
   // === Initialization
   void setCurvesFromEdges(ManifoldSurfaceMesh& mesh);
 
@@ -80,6 +80,9 @@ class NormalCoordinates {
 
   std::array<int, 3> computeBoundaryEdgeSplitDataGeodesic(IntrinsicGeometryInterface& geo, Edge e, double location);
 
+  // Compute the new normal coordinates after splitting edge e
+  // Edge e is split at segment iSeg, where 0 <= iSeg <= edgeCoords[e]
+  std::array<int, 4> computeInteriorEdgeSplitDataCombinatorial(IntrinsicGeometryInterface& geo, Edge e, size_t iSeg);
 
   // Check that these are valid normal coordinates, throw if not
   void validate() const;
@@ -200,7 +203,8 @@ std::vector<SurfacePoint> generateSingleGeodesicGeometry(ManifoldSurfaceMesh& me
 // Barycentric coordinates are 0 at the src of and edge and 1 at the dst
 std::vector<std::pair<SurfacePoint, double>> generateFullSingleGeodesicGeometry(ManifoldSurfaceMesh& mesh,
                                                                                 IntrinsicGeometryInterface& geo,
-                                                                                const NormalCoordinatesCurve& curve);
+                                                                                const NormalCoordinatesCurve& curve,
+                                                                                bool verbose = false);
 
 // Compute the new normal coordinates for an inserted vertex v in face f given
 // that v's barycentric coordinates and the barycentric coordinates of the

include/geometrycentral/surface/signpost_intrinsic_triangulation.h

@@ -49,6 +49,8 @@ class SignpostIntrinsicTriangulation : public IntrinsicTriangulation {
   std::vector<SurfacePoint> traceInputHalfedgeAlongIntrinsic(Halfedge inputHe) override;
   std::vector<SurfacePoint> traceInputHalfedgeAlongIntrinsic(Halfedge inputHe, bool trimEnd);
 
+  bool checkEdgeOriginal(Edge e) const override;
+
   // ======================================================
   // ======== Low-Level Mutators
   // ======================================================
@@ -109,10 +111,6 @@ class SignpostIntrinsicTriangulation : public IntrinsicTriangulation {
   // Scale factor to take Euclidean data to cone data
   double vertexAngleScaling(Vertex v) const;
 
-  // Recover t-values after tracing
-  // Note that really we ought to just report these back from the tracing routine itself, which computes them internally. We don't have a nice API for passing that data around, so this lazily recovers it
-  std::vector<double> recoverTraceTValues(const std::vector<SurfacePoint>& edgeTrace);
-
   // Construct the common subdivision for the current triangulation.
   // WARNING: this signpost implementation does not properly populate some fields
   // of the common subdivision, in particular the `orientation` field for intersections.

include/geometrycentral/surface/surface_point.h

@@ -19,11 +19,11 @@ enum class SurfacePointType { Vertex = 0, Edge, Face };
 struct SurfacePoint {
 
   // === Constructors
-  SurfacePoint();                           // default: yields invalid SurfacePoint with Type::Vertex and null vertex
-  SurfacePoint(Vertex v);                   // at vertex
-  SurfacePoint(Edge e, double tEdge);       // in edge
+  SurfacePoint();                              // default: yields invalid SurfacePoint with Type::Vertex and null vertex
+  SurfacePoint(Vertex v);                      // at vertex
+  SurfacePoint(Edge e, double tEdge);          // in edge
   SurfacePoint(Halfedge he, double tHalfedge); // in edge (flips direction if needed)
-  SurfacePoint(Face f, Vector3 faceCoords); // in face
+  SurfacePoint(Face f, Vector3 faceCoords);    // in face
 
 
   // === Identifying data
@@ -52,6 +52,10 @@ struct SurfacePoint {
   // on or adjacent to the requested face, throws an error.
   inline SurfacePoint inFace(Face f) const;
 
+  // Returns the surface point as an edge point in edge e. If the the SurfacePoint is not on the edge or one of its
+  // adjacent vertices, throws an error
+  inline SurfacePoint inEdge(Edge e) const;
+
   // Return the nearest vertex to this surface point
   inline Vertex nearestVertex() const;
 
@@ -92,4 +96,3 @@ std::string to_string(geometrycentral::surface::SurfacePoint p);
 }
 
 #include "geometrycentral/surface/surface_point.ipp"
-

include/geometrycentral/surface/surface_point.ipp

@@ -147,6 +147,27 @@ inline SurfacePoint SurfacePoint::inFace(Face targetFace) const {
   return *this;
 }
 
+inline SurfacePoint SurfacePoint::inEdge(Edge targetEdge) const {
+  switch (type) {
+  case SurfacePointType::Vertex:
+    if (vertex == targetEdge.halfedge().tailVertex()) {
+      return SurfacePoint(targetEdge, 0);
+    } else if (vertex == targetEdge.halfedge().tipVertex()) {
+      return SurfacePoint(targetEdge, 1);
+    }
+    break;
+  case SurfacePointType::Edge:
+    if (edge == targetEdge) {
+      return *this;
+    }
+    break;
+  default:
+    break;
+  }
+  throw std::logic_error("SurfacePoint " + std::to_string(*this) + " not adjacent to target edge " +
+                         std::to_string(targetEdge));
+  return *this;
+}
 
 inline Vertex SurfacePoint::nearestVertex() const {
 
@@ -292,4 +313,3 @@ inline std::string to_string(geometrycentral::surface::SurfacePoint p) {
   return output.str();
 }
 } // namespace std
-