Signed Heat Method from curve: Each Curve Segment must be contained within a single face ERROR

[Hjamesforsythe]

I am attempting to run a signed heat geodesic distance but failing because of a:

"RuntimeError: Each curve segment must be contained within a single face." Error

The current code looks like this

V, F = mesh_to_VF(mesh)
curves = polylines_to_bary_curves(mesh, polylines)
is_signed = coerce_is_signed(signed if 'signed' in globals() else True, len(curves))
pts_unsigned = sanitize_points_list(source_vertices if 'source_vertices' in globals() else None, V.shape[0])

solver = pp3d.MeshSignedHeatSolver(V, F)
dist = solver.compute_distance(curves, is_signed, pts_unsigned)
dist = np.asarray(dist, dtype=float).ravel().tolist()

I am running this in rhino/grasshopper so the generation of the curves list looks like this:

def _bary_from_meshpoint(mp):
    """
    Extract (face_index, [b0, b1, b2]) from a Rhino MeshPoint.
    """
    if mp is None or mp.FaceIndex < 0:
        raise ValueError("Invalid MeshPoint (no face hit).")

    fi = int(mp.FaceIndex)
    T = list(mp.T)  # MeshPoint.T is double[]; convert to a Python list
    T = T[:3]
    return fi, T

def polylines_to_bary_curves(m, polys):
    """Project each polyline vertex to mesh and encode as (face_index, [b0,b1,b2]) tuples."""
    curves = []
    if not polys:
        return curves

    for obj in polys:
        # Ensure we have an rg.Polyline to iterate
        pl = None
        if isinstance(obj, rg.Polyline):
            pl = obj
        elif isinstance(obj, rg.PolylineCurve):
            pl = obj.ToPolyline()
            
        if pl is None or pl.Count == 0:
            continue

        nodes = []
        for i in range(pl.Count):
            p = pl[i]  # Polyline supports indexer access
            mp = m.ClosestMeshPoint(p, 1e30)
            if mp is None or mp.FaceIndex < 0:
                continue
            fi, b = _bary_from_meshpoint(mp)
            nodes.append((fi, b))
        if nodes:
            curves.append(nodes)
    return curves

does anyone know where this error is coming from and how to address it?

[nmwsharp]

CC @nzfeng

[Hjamesforsythe]

following up looking at the code in geometry central, is it possible to create a list of barycentric face points with these conditions?

void SignedHeatSolver::buildSignedCurveSource(const Curve& curve, Vector<std::complex<double>>& X0) const {

  size_t nNodes = curve.nodes.size();
  for (size_t i = 0; i < nNodes - 1; i++) {
    const SurfacePoint& pA = curve.nodes[i];
    const SurfacePoint& pB = curve.nodes[i + 1];
    Edge commonEdge = sharedEdge(pA, pB);
    if (commonEdge != Edge()) {
      size_t eIdx = geom.edgeIndices[commonEdge];
      std::complex<double> innerProd(0, 1);
      if (pA.vertex == commonEdge.secondVertex()) innerProd.imag(-1);
      double length = lengthOfSegment(pA, pB);
      std::complex<double> contrib = length * innerProd;
      X0[eIdx] += contrib;
      continue;
    }
    Face commonFace = sharedFace(pA, pB);
    if (commonFace == Face()) {
      throw std::logic_error("Each curve segment must be contained within a single face.");
    }
    for (Edge e : commonFace.adjacentEdges()) {
      size_t eIdx = geom.edgeIndices[e];
      std::complex<double> innerProd = projectedNormal(pA, pB, e); // already incorporates length
      X0[eIdx] += innerProd;
    }
  }
}

[nzfeng]

Hi, thanks for your interest! (thanks for cc'ing me, Nick!)

The issue is that the vertices of your polyline, once projected onto the surface mesh, are not guaranteed to form a fully-defined path on the surface mesh:

In particular, adjacent vertices of the polyline could get projected, for example, to faces that are not adjacent in the mesh, leaving an ambiguity of how to connect up the projected points into curves. So the signed heat method code throws an error.

A reasonable thing to do is to simply connect projected points via geodesic paths, perhaps using the EdgeFlipGeodesicSolver described here. In principle, this solver can give you exactly the straightest path between two surface points, though it looks like the solver currently only supports start/end points which are vertices, not general SurfacePoints... is this correct @nmwsharp?

(Might be worth incorporating connecting-points-into-curves-via-geodesic-paths as a feature in the signed heat method code, to be done automatically if the input curves are not fully specified like so...)

[Hjamesforsythe]

hey @nzfeng thanks for getting back to me, I was able to get something to work by sending a long list of curves each defined as two points in the same face. looking something like this:

this is kind of a work around and definitely lead to a more noise than it would lead to with a single fully connected curve.

i have yet to try the geodesic edge flips between consecutive points but as you mentioned when i read the documentation/implementation only works for vertex sources right now so it may not be directly helpful yet

[Hjamesforsythe]

I was able to set up this test case that puts one point in each face that is adjacent to the next point in the next face:

and then checked if each point is adjacent in using face<->face adjacency check:

and still was getting the issue :/

I think that it is coming from this part of the code in geomtrycentral:

 for (size_t i = 0; i < nNodes - 1; i++) {
    const SurfacePoint& pA = curve.nodes[i];
    const SurfacePoint& pB = curve.nodes[i + 1];
...
...
    Face commonFace = sharedFace(pA, pB);
    if (commonFace == Face()) {
      throw std::logic_error("Each curve segment must be contained within a single face.");
    }

in the "SignedHeatSolver::buildSignedCurveSource" function where sharedface(pA,pB) checks:

inline Face sharedFace(const BarycentricVector& u, const BarycentricVector& w) {

  switch (u.type) {
  case BarycentricVectorType::Face: {
    switch (w.type) {
    case BarycentricVectorType::Face:
      if (u.face == w.face) return u.face;
      break;

if they are both faces. Meaning all consecutive face points have to fall on the same face, do you think this is where this issue could be coming from?

[nzfeng]

Sorry, I should have fully specified: the path of a curve along a mesh must include all edge crossings as well: if two surface points are merely in adjacent faces, then there is still an ambiguity of where, exactly, the curve crosses their shared edge, and hence what the resulting two segments of the curve are along the surface -- the geometry of these segments determines the normals of those curve segments, and those normals are used by the signed heat method.

Also, the relevant sharedFace function here is the one for SurfacePoints rather than BarycentricVector (confusing, I know!), whose implementation is here -- and indeed checks if two SurfacePoints lie within a common face.

Appreciate you taking the initiative on these tests!

[Hjamesforsythe]

one of the issues with using the edgeflipsgeodesics -> signedheatmethod pipeline here is that the edgeflipgeodesics returns 3d positions and the signedheatmethod takes in a list of surface points meaning they don't have a direect way of talking to eachother and being used in sequence.

is there a way to return the surface points from the edgeflips mechanism or have the signedheatmethod take 3d positions and do the projection and walking in itself, that would make that workflow accessible

Hi, thanks for your interest! (thanks for cc'ing me, Nick!)

The issue is that the vertices of your polyline, once projected onto the surface mesh, are not guaranteed to form a fully-defined path on the surface mesh:

In particular, adjacent vertices of the polyline could get projected, for example, to faces that are not adjacent in the mesh, leaving an ambiguity of how to connect up the projected points into curves. So the signed heat method code throws an error.

A reasonable thing to do is to simply connect projected points via geodesic paths, perhaps using the EdgeFlipGeodesicSolver described here. In principle, this solver can give you exactly the straightest path between two surface points, though it looks like the solver currently only supports start/end points which are vertices, not general SurfacePoints... is this correct @nmwsharp?

(Might be worth incorporating connecting-points-into-curves-via-geodesic-paths as a feature in the signed heat method code, to be done automatically if the input curves are not fully specified like so...)

[nzfeng]

You're right, it looks like potpourri3d currently doesn't have the option to return geodesic paths as a sequence of SurfacePoints.

But, there's no reason it can't! The geometry-central, the C++ library that potpourri3d contains bindings to, has a pretty complete suite of tools for constructing geodesic paths, which includes returning the path as a sequence of SurfacePoints. It's just that potpourri3d doesn't currently expose a lot of them.

When I get some spare time, I'll make a PR to geometry-central so that the signed heat method pre-processes curves to fill in gaps with geodesic paths, and update potpourri3d. In principle, all you need to do is insert code at line 34 of geometry-central/src/surface/signed_heat_method.cpp so that it looks something like this:

VertexData<double> SignedHeatSolver::computeDistance(const std::vector<Curve>& curves_,
                                                     const std::vector<SurfacePoint>& points,
                                                     const SignedHeatOptions& options) {

  std::vector<Curve> curves = preprocessCurves(curves_);
  Vector<std::complex<double>> Xt = integrateVectorHeatFlow(curves, points, options);
  return integrateVectorField(Xt, curves, points, options);
}

where preprocessCurves is a function that creates an FlipEdgeNetwork that constructs a geodesic path between two SurfacePoints adjacent in a curve if needed (by creating an edge path and then straightening it). I'll try to get to this by next week at the latest...

[Hjamesforsythe]

hey @nzfeng , I come from using geometrycentral very extensively for about 3 years now, and extending/modifying it for my specific use cases in C++ as well. I am very interested in potpourri3d and geometry central continued development / contributing as well!

[nzfeng]

Heads up I've pushed an update to the C++ version in geometry-central here. Upon more thought, I decided to avoid *adding* curves, because in this most generic case, this imposes additional structure that the user might not have intended. Instead, the curve is partitioned in "well-sampled" components ("well-sampled" in the sense I described above). I haven't updated potpourri3d yet.

For your particular use-case -- projecting points onto a mesh -- you'll likely still miss all edge-crossings, so the preprocess step I described above probably won't help much. However, in your particular case, your projected points are at least often projected into adjacent faces...

If you're able to use C++, the following function (adjacentFaces) should check if two SurfacePoints are in at least adjacent faces, and if so, return the common edge between these two faces. Once you have this edge, then you should be able to determine an "optimal" edge-crossing between the two SurfacePoints (what "optimal" means is up to you). Not sure how much this helps...

  // If two SurfacePoints are in a common or adjacent faces, return the faces.
  auto adjacentFaces = [](const SurfacePoint& pA, const SurfacePoint& pB, Edge& commonEdge) -> std::pair<Face, Face> {
    switch (pA.type) {
    case SurfacePointType::Vertex:
      switch (pB.type) {
      case SurfacePointType::Vertex:
        for (Face f : pA.vertex.adjacentFaces()) {
          for (Edge e : f.adjacentEdges()) {
            for (Face g : e.adjacentFaces()) {
              for (Vertex v : g.adjacentVertices()) {
                if (v == pB.vertex) {
                  commonEdge = e;
                  return std::make_pair(f, g);
                }
              }
            }
          }
        }
        break;
      case SurfacePointType::Edge:
        for (Face f : pA.vertex.adjacentFaces()) {
          for (Edge e : f.adjacentEdges()) {
            for (Face g : e.adjacentFaces()) {
              for (Edge h : g.adjacentEdges()) {
                if (h == pB.edge) {
                  commonEdge = e;
                  return std::make_pair(f, g);
                }
              }
            }
          }
        }
        break;
      case SurfacePointType::Face:
        for (Face f : pA.vertex.adjacentFaces()) {
          for (Edge e : f.adjacentEdges()) {
            for (Face g : e.adjacentFaces()) {
              if (g == pB.face) {
                commonEdge = e;
                return std::make_pair(f, g);
              }
            }
          }
        }
        break;
      }
      break;
    case SurfacePointType::Edge:
      switch (pB.type) {
      case SurfacePointType::Vertex:
        for (Face f : pA.edge.adjacentFaces()) {
          for (Edge e : f.adjacentEdges()) {
            for (Face g : e.adjacentFaces()) {
              for (Vertex v : g.adjacentVertices()) {
                if (v == pB.vertex) {
                  commonEdge = e;
                  return std::make_pair(f, g);
                }
              }
            }
          }
        }
        break;
      case SurfacePointType::Edge:
        for (Face f : pA.edge.adjacentFaces()) {
          for (Edge e : f.adjacentEdges()) {
            for (Face g : e.adjacentFaces()) {
              for (Edge h : g.adjacentEdges()) {
                if (h == pB.edge) {
                  commonEdge = e;
                  return std::make_pair(f, g);
                }
              }
            }
          }
        }
        break;
      case SurfacePointType::Face:
        for (Face f : pA.edge.adjacentFaces()) {
          for (Edge e : f.adjacentEdges()) {
            for (Face g : e.adjacentFaces()) {
              if (g == pB.face) {
                commonEdge = e;
                return std::make_pair(f, g);
              }
            }
          }
        }
        break;
      }
      break;
    case SurfacePointType::Face:
      switch (pB.type) {
      case SurfacePointType::Vertex:
        for (Edge e : pA.face.adjacentFaces()) {
          for (Face g : e.adjacentFaces()) {
            for (Vertex v : g.adjacentVertices()) {
              if (v == pB.vertex) {
                commonEdge = e;
                return std::make_pair(pA.face, g);
              }
            }
          }
        }
        break;
      case SurfacePointType::Edge:
        for (Edge e : pA.face.adjacentFaces()) {
          for (Face g : e.adjacentFaces()) {
            for (Edge h : g.adjacentEdges()) {
              if (h == pB.edge) {
                commonEdge = e;
                return std::make_pair(pA.face, g);
              }
            }
          }
        }
        break;
      case SurfacePointType::Face:
        for (Edge e : pA.face.adjacentFaces()) {
          for (Face g : e.adjacentFaces()) {
            if (g == pB.face) {
              commonEdge = e;
              return std::make_pair(pA.face, g);
            }
          }
        }
        break;
      }
      break;
    }
    commonEdge = Edge();
    return std::make_pair(Face(), Face());
  };