Calling FlipEdgeNetwork::delaunayRefine() doesn't seem to preserve the shape of the original mesh

Start by reproducing the reported result with GeodesicPathExtractor::getPathNetworkOnDelaunayMesh() and inspect FlipEdgeNetwork::delaunayRefine(), the intrinsic mesh, and vertexLocations. Compare the refined mesh and interpolated positions with the original image and determine whether the observed edge changes are expected. Done means the behavior is explained and, if it is a bug, a focused regression test and correction are identified.

I'm trying to compute a new mesh from the intrinsic triangulation, but results are not quite as expected. Any clue why some triangles near right edges are being modified like this? Below is the code I'm using to generate a new mesh on which the geodesic path are all surface points on vertices.

PathNetworkSPtr GeodesicPathExtractor::getPathNetworkOnDelaunayMesh(double areaThresh,
                                                                    size_t maxInsertions,
                                                                    double angleBound) const
{
    using namespace geometrycentral::surface;

    //_edgeNetwork->makeDelaunay();
    _edgeNetwork->delaunayRefine(areaThresh, maxInsertions, angleBound);

    auto & tri = _edgeNetwork->tri;
    auto newMeshCopy = tri->intrinsicMesh->copy();

    const auto & vertexData = _mesh->getGeometry().vertexPositions;
    auto * intrinsicMeshPtr = tri->intrinsicMesh.get();
    auto * meshCopyPtr = newMeshCopy.get();

    VertexData<Vector3> newPositions(*newMeshCopy);
    for (auto v : newMeshCopy->vertices()) {
        Vertex vi(intrinsicMeshPtr, v.getIndex());
        auto spOnOrig = tri->vertexLocations[vi];
        newPositions[v] = spOnOrig.interpolate(vertexData);
    }

    const auto vertexIndices = tri->intrinsicMesh->getVertexIndices();

    std::vector<std::vector<SurfacePoint>> newPathNetwork;
    newPathNetwork.reserve(_edgeNetwork->paths.size());
    for (auto & path : _edgeNetwork->paths) {
        const auto heList = path->getHalfedgeList();
        std::vector<SurfacePoint> newSubpath;
        newSubpath.reserve(heList.size() + 1);
        for (auto he : heList) {
            auto vA = Vertex(meshCopyPtr, vertexIndices[he.tailVertex()]);
            auto vB = Vertex(meshCopyPtr, vertexIndices[he.tipVertex()]);
            if (newSubpath.empty() || newSubpath.back() != vA) {
                newSubpath.emplace_back(vA);
            }
            newSubpath.emplace_back(vB);
        }
        newPathNetwork.emplace_back(std::move(newSubpath));
    }

    newMeshCopy->compress();

    auto newGeometry = std::make_unique<VertexPositionGeometry>(*newMeshCopy, newPositions);
    Mesh3DBaseSPtr newMesh = std::make_shared<Mesh3DBase>(std::move(newMeshCopy), std::move(newGeometry));
    return std::make_shared<PathNetwork>(std::move(newPathNetwork), newMesh);
}