Wrong path is generated by FlipEdgeNetwork::iterativeShorten

Start with the provided reproduction and buste.zip, then inspect FlipEdgeNetwork::iterativeShorten and the UV-metric setup in the example. Run the reversed hard-coded path case and compare the generated path with the reported input and UV lengths. Done means identifying why the expected geodesic path is not produced and documenting or correcting the behavior.

I'm calculating the geodesic path of a given (reversed) hard-coded path using UV metric. The path is wrong.

Model:
buste.zip

Code:

#include <algorithm>
#include <cmath>
#include <iostream>
#include <memory>
#include <string>
#include <vector>

#include "geometrycentral/surface/edge_length_geometry.h"
#include "geometrycentral/surface/flip_geodesics.h"
#include "geometrycentral/surface/manifold_surface_mesh.h"
#include "geometrycentral/surface/mesh_graph_algorithms.h"
#include "geometrycentral/surface/meshio.h"
#include "geometrycentral/surface/vertex_position_geometry.h"

namespace GC = geometrycentral;
namespace GCS = geometrycentral::surface;

// Converts a list of contiguous vertex indices into a Halfedge path
std::vector<GCS::Halfedge> build_halfedge_path(
    GCS::ManifoldSurfaceMesh &mesh,
    const std::vector<size_t> &vIndices ) {
    std::vector<GCS::Halfedge> hePath;
    if ( vIndices.size() < 2 )
        return hePath;

    for ( size_t i = 0; i < vIndices.size() - 1; ++i ) {
        GCS::Vertex v1 = mesh.vertex( vIndices[i] );
        GCS::Vertex v2 = mesh.vertex( vIndices[i + 1] );

        GCS::Halfedge targetHe;
        for ( GCS::Halfedge he : v1.outgoingHalfedges() ) {
            if ( he.tipVertex() == v2 ) {
                targetHe = he;
                break;
            }
        }

        if ( !targetHe.getMesh() ) {
            std::cerr << "Error: No direct halfedge found from V" << vIndices[i]
                      << " to V" << vIndices[i + 1] << "\n";
            return {};
        }
        hePath.push_back( targetHe );
    }
    return hePath;
}

// Helper to find candidate faces for a SurfacePoint
std::vector<GCS::Face> getCandidateFaces( const GCS::SurfacePoint &p ) {
    std::vector<GCS::Face> faces;
    if ( p.type == GCS::SurfacePointType::Face ) {
        faces.push_back( p.face );
    } else if ( p.type == GCS::SurfacePointType::Edge ) {
        GCS::Halfedge he = p.edge.halfedge();
        if ( he.isInterior() )
            faces.push_back( he.face() );
        if ( he.twin().isInterior() )
            faces.push_back( he.twin().face() );
    } else if ( p.type == GCS::SurfacePointType::Vertex ) {
        for ( GCS::Face f : p.vertex.adjacentFaces() ) {
            faces.push_back( f );
        }
    }
    return faces;
}

// Helper to find a shared face between two adjacent SurfacePoints
GCS::Face findSharedFace( const GCS::SurfacePoint &p1, const GCS::SurfacePoint &p2 ) {
    auto faces1 = getCandidateFaces( p1 );
    auto faces2 = getCandidateFaces( p2 );
    for ( GCS::Face f1 : faces1 ) {
        for ( GCS::Face f2 : faces2 ) {
            if ( f1 == f2 )
                return f1;
        }
    }
    return GCS::Face();
}

// Helper to get barycentric coordinates of SurfacePoint relative to face f
GC::Vector3 getBarycentricInFace( const GCS::SurfacePoint &p, GCS::Face f ) {
    GCS::Halfedge he0 = f.halfedge();
    GCS::Halfedge he1 = he0.next();
    GCS::Halfedge he2 = he1.next();

    if ( p.type == GCS::SurfacePointType::Vertex ) {
        if ( p.vertex == he0.vertex() )
            return { 1.0, 0.0, 0.0 };
        if ( p.vertex == he1.vertex() )
            return { 0.0, 1.0, 0.0 };
        if ( p.vertex == he2.vertex() )
            return { 0.0, 0.0, 1.0 };
    } else if ( p.type == GCS::SurfacePointType::Edge ) {
        if ( p.edge == he0.edge() ) {
            double t = ( p.edge.halfedge() == he0 ) ? p.tEdge : ( 1.0 - p.tEdge );
            return { 1.0 - t, t, 0.0 };
        }
        if ( p.edge == he1.edge() ) {
            double t = ( p.edge.halfedge() == he1 ) ? p.tEdge : ( 1.0 - p.tEdge );
            return { 0.0, 1.0 - t, t };
        }
        if ( p.edge == he2.edge() ) {
            double t = ( p.edge.halfedge() == he2 ) ? p.tEdge : ( 1.0 - p.tEdge );
            return { t, 0.0, 1.0 - t };
        }
    } else if ( p.type == GCS::SurfacePointType::Face ) {
        if ( p.face == f )
            return p.faceCoords;
    }
    return { 0.0, 0.0, 0.0 };
}

// Seam-safe local face UV interpolation
GC::Vector2 get_uv_in_face( const GCS::SurfacePoint &p, GCS::Face f, const GCS::CornerData<GC::Vector2> &uvs ) {
    GCS::Halfedge he0 = f.halfedge();
    GCS::Halfedge he1 = he0.next();
    GCS::Halfedge he2 = he1.next();

    GC::Vector3 b = getBarycentricInFace( p, f );
    return b.x * uvs[he0.corner()] + b.y * uvs[he1.corner()] + b.z * uvs[he2.corner()];
}

int main() {
    bool reversePath = 1;

    std::string meshFilename = "c:/prj/test_data/relative/_tmp/8/8a/buste.obj";
    std::vector<size_t> pathVerts = { 46961, 46962, 46963, 50001, 49999, 47345, 50000, 47346, 49998, 46987, 46988 };

    if ( reversePath ) {
        std::cout << "[Info] Reversing input path sequence...\n";
        std::reverse( pathVerts.begin(), pathVerts.end() );
    }

    // 1. Load mesh, 3D positions, and corner UVs
    std::cout << "Loading mesh and UV coordinates: " << meshFilename << "...\n";
    std::unique_ptr<GCS::ManifoldSurfaceMesh> mesh;
    std::unique_ptr<GCS::VertexPositionGeometry> geometry;
    std::unique_ptr<GCS::CornerData<GC::Vector2>> uvs;

    std::tie( mesh, geometry, uvs ) = GCS::readParameterizedManifoldSurfaceMesh( meshFilename );

    if ( !uvs ) {
        std::cerr << "Error: OBJ file contains no UV coordinates.\n";
        return 1;
    }

    // 2. Compute 2D UV metric edge lengths for all mesh edges
    GCS::EdgeData<double> uvEdgeLengths( *mesh );
    for ( GCS::Edge e : mesh->edges() ) {
        GCS::Halfedge he = e.halfedge();
        GC::Vector2 uvTail = ( *uvs )[he.corner()];
        GC::Vector2 uvTip = ( *uvs )[he.next().corner()];
        uvEdgeLengths[e] = ( uvTip - uvTail ).norm();
    }

    // Wrap UV metric into EdgeLengthGeometry for FlipEdgeNetwork
    GCS::EdgeLengthGeometry uvGeometry( *mesh, uvEdgeLengths );

    // 3. Build input path
    std::vector<GCS::Halfedge> hePath = build_halfedge_path( *mesh, pathVerts );
    if ( hePath.empty() ) {
        std::cerr << "Failed to construct valid halfedge path.\n";
        return 1;
    }

    // 4. Print Input Path Edge Lengths
    std::cout << "\n=== Input Path Edge Lengths ===\n";
    std::cout << "Seg\tFrom -> To\t\tSurface Length (3D)\tUV Length (2D)\n";
    std::cout << "----------------------------------------------------------------------\n";

    double totalInputSurfaceLength = 0.0;
    double totalInputUVLength = 0.0;

    geometry->requireEdgeLengths();

    for ( size_t i = 0; i < hePath.size(); ++i ) {
        GCS::Halfedge he = hePath[i];
        GCS::Edge e = he.edge();

        double surfaceLen = geometry->edgeLengths[e];
        double uvLen = uvEdgeLengths[e];

        totalInputSurfaceLength += surfaceLen;
        totalInputUVLength += uvLen;

        std::cout << i + 1 << "\tV" << he.tailVertex().getIndex()
                  << " -> V" << he.tipVertex().getIndex()
                  << "\t\t" << surfaceLen
                  << "\t\t" << uvLen << "\n";
    }

    std::cout << "----------------------------------------------------------------------\n";
    std::cout << "Total Input Surface Length (3D): " << totalInputSurfaceLength << "\n";
    std::cout << "Total Input UV Length (2D):      " << totalInputUVLength << "\n\n";

    // 5. Initialize and run geodesic shortening using UV Metric
    std::cout << "Running FlipEdgeNetwork iterativeShorten() using UV Metric...\n";
    std::vector<std::vector<GCS::Halfedge>> gcPaths = { hePath };

    GCS::FlipEdgeNetwork edgeNet( *mesh, uvGeometry, gcPaths );
    edgeNet.posGeom = geometry.get(); // Kept for 3D coordinate evaluation
    edgeNet.addAllWedgesToAngleQueue();

    try {
        edgeNet.iterativeShorten();
        edgeNet.validate();
    } catch ( const std::exception &e ) {
        std::cerr << "GeometryCentral threw exception: " << e.what() << "\n";
        return 1;
    }

    // 6. Inspect geodesic path points and print edge lengths
    std::vector<std::vector<GCS::SurfacePoint>> refinedPaths = edgeNet.getPathPolyline();
    const auto &geodesicPath = refinedPaths[0];

    std::cout << "\n=== Shortened Geodesic Path Points (" << geodesicPath.size() << " points) ===\n";
    for ( size_t i = 0; i < geodesicPath.size(); ++i ) {
        const auto &sp = geodesicPath[i];
        if ( sp.type == GCS::SurfacePointType::Vertex ) {
            std::cout << "  [" << i << "] Vertex V" << sp.vertex.getIndex() << "\n";
        } else if ( sp.type == GCS::SurfacePointType::Edge ) {
            std::cout << "  [" << i << "] Edge E" << sp.edge.getIndex()
                      << " (t=" << sp.tEdge << ")\n";
        } else if ( sp.type == GCS::SurfacePointType::Face ) {
            std::cout << "  [" << i << "] Face F" << sp.face.getIndex()
                      << " (coords=[" << sp.faceCoords.x << ", "
                      << sp.faceCoords.y << ", "
                      << sp.faceCoords.z << "])\n";
        }
    }

    auto formatPointStr = []( const GCS::SurfacePoint &sp ) -> std::string {
        if ( sp.type == GCS::SurfacePointType::Vertex ) {
            return "V" + std::to_string( sp.vertex.getIndex() );
        } else if ( sp.type == GCS::SurfacePointType::Edge ) {
            return "E" + std::to_string( sp.edge.getIndex() ) + "(t=" + std::to_string( sp.tEdge ).substr( 0, 4 ) + ")";
        } else if ( sp.type == GCS::SurfacePointType::Face ) {
            return "F" + std::to_string( sp.face.getIndex() );
        }
        return "Unknown";
    };

    // Seam-safe UV distance calculation
    auto compute_uv_distance = [&]( const GCS::SurfacePoint &p1, const GCS::SurfacePoint &p2 ) -> double {
        GCS::Face f = findSharedFace( p1, p2 );
        if ( !f.getMesh() )
            return 0.0;

        GC::Vector2 uv1 = get_uv_in_face( p1, f, *uvs );
        GC::Vector2 uv2 = get_uv_in_face( p2, f, *uvs );

        return ( uv2 - uv1 ).norm();
    };

    std::cout << "\n=== Shortened Geodesic Path Edge Lengths ===\n";
    std::cout << "Seg\tFrom -> To\t\t\t\tSurface Length (3D)\tUV Length (2D)\n";
    std::cout << "--------------------------------------------------------------------------------\n";

    double totalGeodesicSurfaceLength = 0.0;
    double totalGeodesicUVLength = 0.0;

    for ( size_t i = 0; i < geodesicPath.size() - 1; ++i ) {
        const GCS::SurfacePoint &p1 = geodesicPath[i];
        const GCS::SurfacePoint &p2 = geodesicPath[i + 1];

        // 3D position distance
        GC::Vector3 pos1 = p1.interpolate( geometry->inputVertexPositions );
        GC::Vector3 pos2 = p2.interpolate( geometry->inputVertexPositions );
        double surfaceLen = ( pos2 - pos1 ).norm();

        // Seam-safe local 2D UV distance
        double uvLen = compute_uv_distance( p1, p2 );

        totalGeodesicSurfaceLength += surfaceLen;
        totalGeodesicUVLength += uvLen;

        std::string p1Label = formatPointStr( p1 );
        std::string p2Label = formatPointStr( p2 );

        std::cout << i + 1 << "\t" << p1Label << " -> " << p2Label
                  << "\t\t" << surfaceLen
                  << "\t\t" << uvLen << "\n";
    }

    std::cout << "--------------------------------------------------------------------------------\n";
    std::cout << "Total Geodesic Surface Length (3D): " << totalGeodesicSurfaceLength << "\n";
    std::cout << "Total Geodesic UV Length (2D):      " << totalGeodesicUVLength << "\n\n";

    return 0;
}

Output:

[Info] Reversing input path sequence...
Loading mesh and UV coordinates: c:/prj/test_data/relative/_tmp/8/8a/buste.obj...

=== Input Path Edge Lengths ===
Seg     From -> To              Surface Length (3D)     UV Length (2D)
----------------------------------------------------------------------
1       V46988 -> V46987                0.533525                0.385054
2       V46987 -> V49998                0.5332          0.494299
3       V49998 -> V47346                0.861287                0.596416
4       V47346 -> V50000                0.719734                0.645875
5       V50000 -> V47345                0.632819                0.487869
6       V47345 -> V49999                0.739344                0.569489
7       V49999 -> V50001                0.863134                0.64117
8       V50001 -> V46963                1.25565         0.808103
9       V46963 -> V46962                0.995753                0.519554
10      V46962 -> V46961                1.04801         0.498308
----------------------------------------------------------------------
Total Input Surface Length (3D): 8.18246
Total Input UV Length (2D):      5.64614

Running FlipEdgeNetwork iterativeShorten() using UV Metric...

=== Shortened Geodesic Path Points (12 points) ===
  [0] Vertex V46988
  [1] Edge E141860 (t=2.94659e-10)
  [2] Edge E141862 (t=1)
  [3] Edge E142989 (t=0.93919)
  [4] Edge E142991 (t=0.944456)
  [5] Edge E144087 (t=0.825262)
  [6] Edge E144089 (t=0.581222)
  [7] Edge E145272 (t=0.244087)
  [8] Edge E146557 (t=0.172933)
  [9] Edge E146558 (t=0.351293)
  [10] Edge E146560 (t=0.680193)
  [11] Face F96859 (coords=[0.66761, 0.0393139, 0.293076])

=== Shortened Geodesic Path Edge Lengths ===
Seg     From -> To                              Surface Length (3D)     UV Length (2D)
--------------------------------------------------------------------------------
1       V46988 -> E141860(t=0.00)               0.342061                0.493417
2       E141860(t=0.00) -> E141862(t=1.00)              2.91937e-10             3.63525e-10
3       E141862(t=1.00) -> E142989(t=0.93)              1.09776         1.36326
4       E142989(t=0.93) -> E142991(t=0.94)              0.0238264               0.0318704
5       E142991(t=0.94) -> E144087(t=0.82)              0.481669                0.619044
6       E144087(t=0.82) -> E144089(t=0.58)              0.355933                0.442551
7       E144089(t=0.58) -> E145272(t=0.24)              0.350507                0.444904
8       E145272(t=0.24) -> E146557(t=0.17)              0.22725         0.263041
9       E146557(t=0.17) -> E146558(t=0.35)              0.422288                0.451218
10      E146558(t=0.35) -> E146560(t=0.68)              0.671649                0.814278
11      E146560(t=0.68) -> F96859               0.388592                0.387229
--------------------------------------------------------------------------------
Total Geodesic Surface Length (3D): 4.36154
Total Geodesic UV Length (2D):      5.31081

The path is completely in the wrong direction, and it ends up in a face.

The nonreversed path (setting reversePath = 0) is fine:

Loading mesh and UV coordinates: c:/prj/test_data/relative/_tmp/8/8a/buste.obj...

=== Input Path Edge Lengths ===
Seg     From -> To              Surface Length (3D)     UV Length (2D)
----------------------------------------------------------------------
1       V46961 -> V46962                1.04801         0.498308
2       V46962 -> V46963                0.995753                0.519554
3       V46963 -> V50001                1.25565         0.808103
4       V50001 -> V49999                0.863134                0.64117
5       V49999 -> V47345                0.739344                0.569489
6       V47345 -> V50000                0.632819                0.487869
7       V50000 -> V47346                0.719734                0.645875
8       V47346 -> V49998                0.861287                0.596416
9       V49998 -> V46987                0.5332          0.494299
10      V46987 -> V46988                0.533525                0.385054
----------------------------------------------------------------------
Total Input Surface Length (3D): 8.18246
Total Input UV Length (2D):      5.64614

Running FlipEdgeNetwork iterativeShorten() using UV Metric...

=== Shortened Geodesic Path Points (18 points) ===
  [0] Vertex V46961
  [1] Edge E141748 (t=0.713315)
  [2] Edge E141750 (t=0.788974)
  [3] Edge E141752 (t=0.760111)
  [4] Edge E142916 (t=0.407001)
  [5] Edge E142915 (t=0.352846)
  [6] Edge E151076 (t=0.669353)
  [7] Edge E151077 (t=0.386777)
  [8] Edge E142948 (t=0.670536)
  [9] Edge E142949 (t=0.112329)
  [10] Edge E142952 (t=0.161902)
  [11] Edge E142954 (t=0.0361707)
  [12] Edge E144061 (t=0.0236971)
  [13] Edge E144063 (t=0.111696)
  [14] Edge E151072 (t=0.606278)
  [15] Edge E141818 (t=0.920299)
  [16] Edge E141819 (t=0.099154)
  [17] Vertex V46988

=== Shortened Geodesic Path Edge Lengths ===
Seg     From -> To                              Surface Length (3D)     UV Length (2D)
--------------------------------------------------------------------------------
1       V46961 -> E141748(t=0.71)               1.17098         0.39314
2       E141748(t=0.71) -> E141750(t=0.78)              0.751442                0.413967
3       E141750(t=0.78) -> E141752(t=0.76)              0.21595         0.136929
4       E141752(t=0.76) -> E142916(t=0.40)              0.413084                0.280327
5       E142916(t=0.40) -> E142915(t=0.35)              0.57216         0.357482
6       E142915(t=0.35) -> E151076(t=0.66)              0.585           0.432296
7       E151076(t=0.66) -> E151077(t=0.38)              0.393859                0.303475
8       E151077(t=0.38) -> E142948(t=0.67)              0.489059                0.374379
9       E142948(t=0.67) -> E142949(t=0.11)              0.279486                0.21315
10      E142949(t=0.11) -> E142952(t=0.16)              0.549951                0.428621
11      E142952(t=0.16) -> E142954(t=0.03)              0.734039                0.579586
12      E142954(t=0.03) -> E144061(t=0.02)              0.0424238               0.0369322
13      E144061(t=0.02) -> E144063(t=0.11)              0.114232                0.091589
14      E144063(t=0.11) -> E151072(t=0.60)              0.773209                0.612231
15      E151072(t=0.60) -> E141818(t=0.92)              0.457032                0.271668
16      E141818(t=0.92) -> E141819(t=0.09)              0.0809719               0.076756
17      E141819(t=0.09) -> V46988               0.488703                0.308287
--------------------------------------------------------------------------------
Total Geodesic Surface Length (3D): 8.11159
Total Geodesic UV Length (2D):      5.31081

Surprsingly, they both have the same length.