Logarithmic map and geodesic tracing for general polygon meshes

Start by reading SurfacePoint, VectorHeatSolver::computeLogMap, and traceGeodesic to identify their barycentric and triangular-mesh assumptions. The issue is complete only when logarithmic maps and geodesic tracing support the listed general polygon cases, including boundaries, non-adjacent faces, shared vertices, multiple shared edges, and concave unfoldings.

My use case

Given a manifold polygon mesh

• possibly with boundary
• with no limit on the number of edges in its faces
• plus two points projected onto the mesh anywhere (not necessarily on vertices or edges)

trace a geodesic path between the two points. If it is not the shortest, at least be straightest (Polthier and Schmies [1998]).

I will use such an algorithm to cut up a surface with profile curves (i.e. curves on a surface). In practice, profile curves are polylines with vertices projected onto the polygon mesh. The idea is that these curve polyline segments become new edges on existing faces, splitting them, while preserving the shape of the original mesh. The curve is discretized finely enough such that most polyline segments have both endpoints inside one face, and these are of no concern. The exceptions are the focus of this issue, when segment endpoints do not reside in one face.

Mock up usage (made with Blender)

Polygon mesh is given (quad mesh in this example), plus the highlighted points. The highlighted points are endpoints of the geodesic trace query, and the query results are the edges drawn across the faces, represented as ordered edge intersection points.

My current but insufficient workarounds

In the simple case that the endpoints of a segment reside on different faces that share exactly one edge (sometimes, see the last failure case below), I work around the problem by calculating where the curve segment would intersect the shared edge if the adjacent faces were unfolded. The intersection point is used to subdivide both the curve segment and the mesh edge so that the curve conforms to the surface. This workaround calculates a straightest geodesic (see first mockup image), but it does not cover all cases. Specifically, it cannot handle the following cases of segment endpoint locations that I have encountered:

• If they end up on non-adjacent polygons
  • This can happen if the segment "skips over" a small gap in the mesh.
  • I do not have a workaround for this case.
• If they end up on adjacent polygons that share a vertex
  • One possible workaround is to connect the endpoints through the mesh vertex, but the path will not be a shortest or straightest geodesic in general.
  • The second mock up image demonstrates this case
• If they end up on adjacent polygons that share more than one edge
  • Which of the shared edges should I unfold along to calculate the intersection?
  • This case can happen when two kite-shaped quads are glued together back-to-back, among other things.
  • My workaround is to just pick "the first one" arbitrarily, but it frequently falls through to the next case.
• If they end up on adjacent polygons that share only one edge, but the unfolded shape is concave
  • Unfolding to calculate the segment-edge intersection fails because the "intersection" would be found beyond one of the ends the mesh edge, outside of both faces. This is in fact the definition of concave.
  • My workaround is to just clamp the point to the mesh edge. Obviously this does not yield a geodesic because it turns into the second case.

How can I help?

If no one is working on such a solution, how can I do it myself? Off the top of my head, I recognize several changes that must be made to implement the algorithm:

• If I am to use generalized barycentric coordinates to represent a point on a face, the SurfacePoint type needs to support varying numbers of components for a face point, depending on the face.
  • I will also need to implement functions that map between generalized barycentric coordinates and 3D Euclidean coordinates.
  • I do not know how to actually do this. I came across a paper by Mark Meyer, Haeyoung Lee, Alan Barr, Mathieu Desbrun [2002], but according to the authors this approach is only applicable to convex polygons, not concave polygons.
• The Vector Heat Method and the VectorHeatSolver::computeLogMap function need to be updated to support general polygonal meshes.
  • Apparently this was already done in the Vector Heat Method paper, in section 6.2.1, but it was not released as part of Geometry Central.
  • The computeLogMap function calculates the log map for a point in a face by taking the barycentric weighted average of the log maps from vertices adjacent to the face. I will need to modify it to use generalized barycentric coordinates or some other interpolation scheme.
  • I am aware of another polygonal discretization of the Laplace-Beltrami operator and other differential operators that has since been published after the Vector Heat Method paper, from Fernando de Goes, Andrew Butts, and Mathieu Desbrun. [2020], but it describes operators on faces and face corners rather than on mesh vertices. I do not know how to properly convert the operators on faces to vertices.
• The traceGeodesic function uses barycentric coordinates for geodesic tracing in triangles, so it needs to be modified to use some other tangent space parameterization to perform the trace in general polygons.
  • I do not know how to do this.

Holdover options

Since this problem is only part of a larger project, I have to choose some other workaround that does not fail in the meantime. I am aware that I can simply triangulate the given mesh, perhaps with constrained Delaunay triangulation, then use the existing computeLogMap and traceGeodesic functions. However, as noted by Fernando de Goes, Andrew Butts, and Mathieu Desbrun. [2020], Figure 2, triangulating a mesh may lead to inaccurate parameterizations, and I would like to avoid that where possible. I am open to suggestions.