SurfacePoint::inEdge Vertex-typed returns reflected tEdge (master); int-tri-updates has a fix

[designbynumbers]

SurfacePoint::inEdge returns reflected tEdge for Vertex-typed sources

I hit this while building a research codebase on top of gc.
Before opening anything I wanted to share what I found and ask how
you'd prefer to handle it, since you've already addressed the same
bug on the int-tri-updates branch (see below) and master still
exhibits it.

The bug

SurfacePoint::inEdge for a Vertex-typed source returns the
reflected edge parameter
(include/geometrycentral/surface/surface_point.ipp:151-157 on
master):

case SurfacePointType::Vertex:
  if (vertex == targetEdge.halfedge().tailVertex()) {
    return SurfacePoint(targetEdge, 1);     // -> reads as tip, should be tail
  } else if (vertex == targetEdge.halfedge().tipVertex()) {
    return SurfacePoint(targetEdge, 0);     // -> reads as tail, should be tip
  }

The tEdge convention used by every other branch of this file is the
opposite. Four witnesses inside surface_point.ipp itself:

Line	Code	Implies
12	tEdge(he == he.edge().halfedge() ? tHalfedge : (1 - tHalfedge))	walking from he.tailVertex() at tHalfedge=0 means tEdge=0 is at tail of edge.halfedge()
67-68	inFace-from-Edge: (1-tEdge, tEdge, 0) when he == targetHe	tEdge=0 puts all weight on inFace.halfedge().vertex() (= edge tail)
180-181	nearestVertex: if (tEdge < .5) return edge.halfedge().vertex();	tEdge < 0.5 rounds toward the tail
211	interpolate: (1-tEdge) * valTail + tEdge * valTip	tEdge=0 ⇒ tail, tEdge=1 ⇒ tip

So the convention seems pretty clear: tEdge=0 at the tail of
edge.halfedge(), tEdge=1 at the tip. Only the Vertex branch of
inEdge disagrees.

IntegerCoordinatesIntrinsicTriangulation::computeEdgeSplitData in
the normalCoordinates < 0 branch (shared / "Case 0a" edges) calls
tail.inEdge(inputEdge).tEdge and tip.inEdge(inputEdge).tEdge and
linearly interpolates between them. Because both endpoints come back
reflected, the new vertex lands at parameter (1 - t) along the
input edge instead of t. Every insertVertex(SurfacePoint(edge, t))
call on a shared intrinsic edge — and every delaunayRefine edge
split as well — places its Steiner at the wrong location, drifting
worst near the endpoints and least near t = 0.5.

While diagnosing this I checked the int-tri-updates branch and
found you'd already corrected it as part of a larger refactor: the
public SurfacePoint::inEdge is removed from surface_point.{h,ipp}
entirely, and the same logic appears as a private lambda inside
IntegerCoordinatesIntrinsicTriangulation::splitBoundaryEdge (around
src/surface/integer_coordinates_intrinsic_triangulation.cpp:1256)
and again inside another mutator (~1442), both with the correct
tailVertex → tEdge=0, tipVertex → tEdge=1 convention.

It looks like that branch never got merged back to master (last
commit 2021-09-21, master last touched it earlier). Anyone building
against current master who calls inEdge — or who calls
insertVertex on a shared intrinsic edge — hits this.

Here's a MWE which shows that the bug still happens on master:

Self-contained single-.cpp reproducer (just interpolate and
inEdge on a one-triangle mesh; no IT machinery needed):

gc_inedge_bug_mwe.cpp (~100 lines)

// Minimum working example: SurfacePoint::inEdge returns reflected tEdge
// for Vertex-typed sources. See the accompanying PR description for
// the convention argument inside surface_point.ipp itself.
//
// Build (link against an existing geometry-central build):
//
//   c++ -std=c++17 -O2 \
//       -I <gc-source>/include \
//       -I <eigen-include-dir> \
//       gc_inedge_bug_mwe.cpp \
//       <gc-build>/lib/libgeometry-central.a \
//       -o gc-mwe
//   ./gc-mwe
//
// Expected output on a vanilla (unpatched) gc:
//
//   Edge-typed (e01, t=0)  interpolates to (0,0,0); expected eTail position (0,0,0)
//   inEdge(Vertex at eTail) -> (e01, tEdge=1) -> 3D (1,0,0); expected eTail position (0,0,0)
//     ** BUG: inEdge(Vertex@tail) returned tEdge=1 (should be 0)
//   inEdge(Vertex at eTip)  -> (e01, tEdge=0) -> 3D (0,0,0); expected eTip position (1,0,0)
//     ** BUG: inEdge(Vertex@tip) returned tEdge=0 (should be 1)
//
// With the two-line fix applied to surface_point.ipp's inEdge
// Vertex-typed branch, both round-trips print `OK: round-trip matches.`

#include <cassert>
#include <iostream>
#include <vector>

#include "geometrycentral/surface/manifold_surface_mesh.h"
#include "geometrycentral/surface/surface_point.h"
#include "geometrycentral/surface/vertex_position_geometry.h"
#include "geometrycentral/utilities/vector3.h"

namespace gcs = geometrycentral::surface;
using geometrycentral::Vector3;

int main() {
    // Single-triangle mesh. Any fixture with at least one edge works;
    // the bug is convention-level, not topology-dependent.
    std::vector<std::vector<size_t>> tris = {{0, 1, 2}};
    auto mesh = std::make_unique<gcs::ManifoldSurfaceMesh>(tris);

    gcs::VertexData<Vector3> pos(*mesh);
    auto vIt = mesh->vertices().begin();
    gcs::Vertex v0 = *vIt; ++vIt;
    gcs::Vertex v1 = *vIt; ++vIt;
    gcs::Vertex v2 = *vIt;
    pos[v0] = Vector3{0.0, 0.0, 0.0};
    pos[v1] = Vector3{1.0, 0.0, 0.0};
    pos[v2] = Vector3{0.0, 1.0, 0.0};

    gcs::VertexPositionGeometry geom(*mesh, pos);

    // Find the edge between v0 and v1.
    gcs::Edge e01;
    for (gcs::Halfedge he : v0.outgoingHalfedges()) {
        if (he.tipVertex() == v1) { e01 = he.edge(); break; }
    }
    assert(e01 != gcs::Edge());

    gcs::Vertex eTail = e01.halfedge().tailVertex();
    gcs::Vertex eTip  = e01.halfedge().tipVertex();

    // Sanity check: the canonical t=0 Edge-typed point interpolates to
    // eTail's position. (This is the tEdge convention surface_point.ipp
    // uses everywhere else; see issue text for the four witnesses.)
    {
        gcs::SurfacePoint sp(e01, 0.0);
        Vector3 p = sp.interpolate(geom.vertexPositions);
        std::cout << "Edge-typed (e01, t=0)  interpolates to ("
                  << p.x << "," << p.y << "," << p.z
                  << "); expected eTail position ("
                  << pos[eTail].x << "," << pos[eTail].y << "," << pos[eTail].z
                  << ")\n";
        assert((p - pos[eTail]).norm() < 1e-12);
    }

    // Bug demonstration: Vertex-typed -> inEdge -> interpolate should
    // round-trip the vertex's position. It does not.
    {
        gcs::SurfacePoint vSP(eTail);
        gcs::SurfacePoint eSP = vSP.inEdge(e01);
        Vector3 p = eSP.interpolate(geom.vertexPositions);
        std::cout << "inEdge(Vertex at eTail) -> (e01, tEdge=" << eSP.tEdge
                  << ") -> 3D (" << p.x << "," << p.y << "," << p.z
                  << "); expected eTail position ("
                  << pos[eTail].x << "," << pos[eTail].y << "," << pos[eTail].z
                  << ")\n";
        if ((p - pos[eTail]).norm() > 1e-12) {
            std::cout << "  ** BUG: inEdge(Vertex@tail) returned tEdge="
                      << eSP.tEdge << " (should be 0)\n";
        } else {
            std::cout << "  OK: round-trip matches.\n";
        }
    }
    {
        gcs::SurfacePoint vSP(eTip);
        gcs::SurfacePoint eSP = vSP.inEdge(e01);
        Vector3 p = eSP.interpolate(geom.vertexPositions);
        std::cout << "inEdge(Vertex at eTip)  -> (e01, tEdge=" << eSP.tEdge
                  << ") -> 3D (" << p.x << "," << p.y << "," << p.z
                  << "); expected eTip position ("
                  << pos[eTip].x << "," << pos[eTip].y << "," << pos[eTip].z
                  << ")\n";
        if ((p - pos[eTip]).norm() > 1e-12) {
            std::cout << "  ** BUG: inEdge(Vertex@tip) returned tEdge="
                      << eSP.tEdge << " (should be 1)\n";
        } else {
            std::cout << "  OK: round-trip matches.\n";
        }
    }

    return 0;
}

Output against current master:

Edge-typed (e01, t=0)  interpolates to (0,0,0); expected eTail position (0,0,0)
inEdge(Vertex at eTail) -> (e01, tEdge=1) -> 3D (1,0,0); expected eTail position (0,0,0)
  ** BUG: inEdge(Vertex@tail) returned tEdge=1 (should be 0)
inEdge(Vertex at eTip)  -> (e01, tEdge=0) -> 3D (0,0,0); expected eTip position (1,0,0)
  ** BUG: inEdge(Vertex@tip) returned tEdge=0 (should be 1)

Can I help somehow? For instance:

1. Minimal master-side PR — a two-line correction to the existing
   inEdge Vertex branch (swap the 0 and 1 returns). Smallest
   possible change for people stuck on master.

2. Help land the last bit of work on int-tri-updates — your
   branch already has the architecturally cleaner fix (remove inEdge
   from public API, private lambdas where it's actually used). If the
   last commits on that branch after it was merged into master are
   just forgotten, I could rebase it onto master, work through any
   conflicts, and put it up as a clean PR.

If you've got a third preference (it's intentional, leave it alone,
known bug, etc.) please just say so-- I have patched my local copy of
gc for now, and I could design around this if needed.

Thanks for gc — it's really fun!

— Jason Cantarella & Claude Opus 4.7