Logarithmic map and geodesic tracing for general polygon meshes

[xprssn2]

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.

[xprssn2]

I see that #195 adds support for polygonal DEC operators and makes the vector heat solver work on polygon meshes. Thanks @nzfeng! That's great because a log map will provide a way (or is at least part of one) to trace geodesic paths across the polygon mesh surface that will connect any two points, with the following procedure:

• Start with the two points on the polygon mesh (may be inside a face, not just on vertices as stated in the original comment)
• Generate a log map using one of the points as the source
  • The triangle mesh vector heat solver does this with barycentric coordinates if the source point is not on a vertex. Will I need to use substitute it with mean value coordinates?
• Sample the log map at the second point, this gives the geodesic distance and tracing direction
• Build a supporting data structure that marks edges that contain the sampled direction within (e.g. if one endpoint of the edge has an angle less than the sample and the other endpoint greater than the sample), and record the inverse-lerp factor for that sampled direction
• Sort the marked edges by geodesic distance
• Then finally, in ascending order of geodesic distance, evaluate the inverse-lerp factor on each edge to get a point within the edge, and chain it one after the next until going through the entire list of marked edges. The result is the geodesic path traced between a source and an end point on a polygon mesh.
  

However

#195 does not (yet?) include an implementation of log map computation for polygon meshes. I have been trying to utilize the polygonal vector heat solver to implement the log map algorithm from the paper by myself, but got stuck. I am not actually trained in DDG so I am not really sure whether what I am trying actually helps. I am more programmer than mathematician.

Here's what I have done so far. (The source point I used in this screenshot is under the wrist, where a watch strap would usually touch.) Transporting a horizontal direction ("H_0" in the Vector Heat Method paper) works flawlessly. This is shown as yellow arrows in my screenshot.

The next step to compute the initial conditions for the radial field ("R" in the original paper, appendix A) is where I got stuck. I do not know how to generalize the hat functions for triangles to polygons. My personal uneducated guess is that the hat function for triangles came from barycentric coordinates (at least the formula for it equals that of the barycentric weight for a triangle vertex). If that is the case, is a generalized barycentric coordinate weight the answer? But if it is, then is there still a closed form solution for the integral of the vertex outward ball? Both Wachspress coordinates and mean value Coordinates are nonlinear and have sums and products and even square roots in the denominators. Will I need a quadrature to calculate the contributions to vertices for the initial condition?

I have also considered another option. Would it be possible to use a gradient operator on the scalar geodesic distance map obtained by scalar heat diffusion, and skip the whole radial field diffusion step by directly calculating the distance gradient? de Goes et al defined a gradient operator for the polygon faces (and I used it above, that's shown as purple arrows), is it possible to come up with one for vertices as well? If yes, how do I do it?

In another attempt, I attempted to simply go through each vertex and add the differences in geodesic distance between its neighbors multiplied by unit vectors pointing from the neighbor back to the vertex itself. This naive version of a "vertex gradient operator" is not correct and does not work, not even crudely.

Any ideas @nmwsharp?