Skip to content

Flip Geodesics

This algorithm takes as input a path (or loop/network of paths) along the edges of a triangle mesh, and as output straightens that path to be a geodesic (i.e. a straight line, or equivalently a locally-shortest path along a surface). The procedure runs in milliseconds, is quite robust, comes with a strong guarantee that no new crossings will be created in the path, and as an added benefit also generates a triangulation on the surface which conforms to the geodesic. Additionally, it even enables the construction of Bézier curves on a surface.

shorten path to geodesic

See this repo for a demo application which wraps this implementation as a standalone application with a GUI.

This is an implementation of the algorithm described in You Can Find Geodesic Paths in Triangle Meshes by Just Flipping Edges.

#include "geometrycentral/surface/flip_geodesics.h"

Example:

#include "geometrycentral/surface/flip_geodesics.h"
#include "geometrycentral/surface/meshio.h"

using namespace geometrycentral;
using namespace geometrycentral::surface;

// Load a mesh
std::unique_ptr<ManifoldSurfaceMesh> mesh;
std::unique_ptr<VertexPositionGeometry> geometry;
std::tie(mesh, geometry) = readManifoldSurfaceMesh("my_mesh.obj");

// Create a path network as a Dijkstra path between endpoints
std::unique_ptr<FlipEdgeNetwork> edgeNetwork;
Vertex vStart = mesh->vertex(0);
Vertex vEnd = mesh->vertex(7);
edgeNetwork = FlipEdgeNetwork::constructFromDijkstraPath(*mesh, *geometry, vStart, vEnd);

// Make the path a geodesic
edgeNetwork->iterativeShorten();

// Extract the result as a polyline along the surface
edgeNetwork->posGeom = geometry.get();
std::vector<std::vector<Vector3>> polyline = edgeNetwork->getPathPolyline3D();

Constructing a path network

The class stateful class FlipEdgeNetwork encapsulates all functionality in the algorithm; all of the methods below are members functions. This class wraps a SignpostIntrinsicTriangulation (as member FlipEdgeNetwork::tri), and also manages a path network along the edges of that triangulation.

The most direct construct takes as input a sequence of halfedges along the surface which comprise the path.

FlipEdgeNetwork::FlipEdgeNetwork(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& inputGeom, const std::vector<std::vector<Halfedge>>& paths, VertexData<bool> extraMarkedVerts = VertexData<bool>())

Construct a network of paths along the surface.

  • mesh: the input manifold triangle mesh
  • inputGeom: the input geometry (as always, a VertexPositionGeometry is valid input)
  • paths: a collection of paths along the surface, each given as an ordered sequence of halfedges
  • extraMarkedVerts: if given, boolean vector of vertices where the path is pinned, and should not be straightened

Additionally, several factory methods are provided which assist in constructing a path network from common input data.

static std::unique_ptr<FlipEdgeNetwork> FlipEdgeNetwork::constructFromDijkstraPath(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, Vertex startVert, Vertex endVert)

Construct a network consisting of a single path, generated by running Dijkstra’s algorithm between two point.

  • mesh: the input manifold triangle mesh
  • inputGeom: the input geometry (as always, a VertexPositionGeometry is valid input)
  • startVert, endVert: the source and target vertices from which to initialize a Dijkstra path along edges
static std::unique_ptr<FlipEdgeNetwork> FlipEdgeNetwork::constructFromPiecewiseDijkstraPath(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, std::vector<Vertex> points, bool closed = false, bool markInterior = false)

Construct a network consisting of a single path, generated by running Dijkstra’s algorithm between consecutive pairs of vertices.

  • mesh: the input manifold triangle mesh
  • inputGeom: the input geometry (as always, a VertexPositionGeometry is valid input)
  • points: the not-necessarily-adjacent vertices from which to initialize the Dijkstra path, by running Dijkstra’s algorithm between each i’th and (i+1)’th points.
  • closed: if true, the first point will be connected to the last to create a closed loop
  • markInterior: if true, the path will be pinned at the specified points and not shortened there
static std::unique_ptr<FlipEdgeNetwork> FlipEdgeNetwork::constructFromEdgeSet(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, const EdgeData<bool>& inPath, const VertexData<bool>& extraMarkedVertices)

Construct a network from a marked set of edges on the mesh. Endpoints, loops, etc are heuristically inferred (e.g. if a path ends at the same vertex it starts at, it is assumed to be a loop).

  • mesh: the input manifold triangle mesh
  • inputGeom: the input geometry (as always, a VertexPositionGeometry is valid input)
  • inPath: a boolean array indicating edges which are in the path
  • extraMarkedVerts: if given, boolean vector of vertices where the path is pinned, and should not be straightened
FlipEdgeNetwork::reinitializePath(const std::vector<std::vector<Halfedge>>& paths)

Re-initialize an existing FlipEdgeNetwork object to hold the paths in paths, just like in the constructor. Useful when used in conjunction with rewind().

  • paths: a collection of paths along the surface, each given as an ordered sequence of halfedges

Manipulating paths

The key routine is FlipEdgeNetwork::iterativeShorten(), which shortens the path until it is a geodesic by flipping edges in the underlying intrinsic triangulation.

void FlipEdgeNetwork::iterativeShorten(size_t maxIterations = INVALID_IND, double maxRelativeLengthDecrease = 0.)

Shorten the network to be geodesic.

Additional termination conditions can be set to potentially stop the process before the path is a geodesic:

  • maxIterations if set to something other than INVALID_IND, limits the number of iterations
  • maxRelativeLengthDecrease if set to something other than 0, limits the maximum decrease in length for the network. E.g. 0.5 would mean the resulting network is at at 0.5 * L, where L is the initial length.

Importantly, after calling iterativeShorten(), the underlying intrinsic triangulation FlipEdgeNetwork::tri is a triangulation which conforms to the generated geodesics. This triangulation can then be used for subsequent computation, such as solving a PDE which has geodesics as boundary conditions.

Shortening many paths

After shortening a path to a geodesic, the underlying triangulation has been mutated by edge flips. If you want to shorten another path, you can call rewind() to reset the data structure and prepare it for the next path.

For instance, if you wanted to find many point-to-point geodesic paths, you might write something like:

std::unique_ptr<FlipEdgeNetwork> flipNetwork(new FlipEdgeNetwork(*mesh, *geometry, {}));
flipNetwork->supportRewinding = true;
flipNetwork->posGeom = geometry.get();

for(/* each path to compute */) {

  Vertex startVert; // populate these somehow
  Vertex endVert; 

  // Get an initial dijkstra path
  // (surface/mesh_graph_algorithms.h)
  std::vector<Halfedge> dijkstraPath = shortestEdgePath(*geom, startVert, endVert);

  // Reinitialize the ede network to contain this path
  flipNetwork->reinitializePath({dijkstraPath});

  // Straighten the path to geodesic
  flipNetwork->iterativeShorten();

  // Extract the path and store it in the vector
  std::vector<Vector3> path3D = flipNetwork->getPathPolyline3D().front();

  // Be kind, rewind
  flipNetwork->rewind();
}

Geodesic Bézier curves

Once we can shorten curves to geodesics, we can use a de Casteljau-style subdivision scheme to construct geodesic Bézier curves along a mesh, using a procedure described in “Modeling on triangulations with geodesic curves” by Morera et al. (2008).

void FlipEdgeNetwork::bezierSubdivide(size_t nRounds)

Construct a Bézier curve from an input control polygon path via iterative subdivision.

  • nRounds specifies the number of subdivision rounds to perform; the curve converges to an exact geodesics as nRounds increases. Setting nRounds=3 is often a good initial choice.

Note that the control points of the curve should be marked vertices, as described in the construction section. Additionally, this function should be applied to a network consisting of a single path, not distinct paths between multiple control points. Calling the factory method FlipEdgeNetwork::constructFromPiecewiseDijkstraPath(mesh, geom, /* vector of constrol vertices */, false, true) above will construct a path ready for Bézier subdivision.

The algorithm used here is essentially the one described in “Modeling on triangulations with geodesic curves” by Morera et al. (2008), but using flip-based geodesics to straighten paths.

A geodesic Bézier curve generated with this method.

Improving the triangulation

As noted above, we can additionally generate a triangulation which conforms to geodesic curves, but this triangulation will often contain extremely skinny triangles. Performing Delaunay refinement improves the triangle quality for subsequent computation.

void FlipEdgeNetwork::delaunayRefine(double areaThresh = std::numeric_limits<double>::infinity(), size_t maxInsertions = INVALID_IND, double angleBound = 25.)

Perform intrinsic Delaunay refinement on the underlying triangulation to improve its numerical quality while preserving the geodesics.

Additional arguments are as in SignpostIntrinsicTriangulation::delaunayRefine().

A geodesic path is introduced, then the triangulation is refined to remove skinny triangles.

Output and utilities

After running iterativeShorten(), the resulting geodesic paths are represented as a collection edges in an updated intrinsic triangulation. The most common operation is to extract the resulting path network as a polyline along the surface. Note that these routines return a double list; the outer list represents the paths/loops in the network, while the inner list are the points in each path/loop.

One option is to get the output as a sequence of surface points, where each is a vertex or a point along some edge in the mesh.

std::vector<std::vector<SurfacePoint>> FlipEdgeNetwork::getPathPolyline()

Get the path network as a sequence of surface points along the input mesh.

std::vector<std::vector<SurfacePoint>> FlipEdgeNetwork::getAllEdgePolyline()

Get all edges in the underlying triangulation as a sequence of surface points along the input mesh.

A simpler approach is to simply extract the path as a sequence of points in 3D space.

To use this approach, you must first set your 3D mesh geometry as the posGeom member. This extra step is necessary because in general, the FlipEdgeNetwork can be applied to meshes which have only edge lengths, but extracting 3D points requires your mesh to have 3D positions. (see the example at the top.)

std::vector<std::vector<Vector3>> FlipEdgeNetwork::getPathPolyline3D()

Get the path network as a sequence of points in 3D space.

To use this approach, you must first set your 3D mesh geometry as the posGeom member.

std::vector<std::vector<Vector3>> FlipEdgeNetwork::getAllEdgePolyline3D()

Get all edges in the underlying intrinsic triangulation as a sequence of points in 3D space.

To use this approach, you must first set your 3D mesh geometry as the posGeom member.

Additionally, a collection of useful utility functions are provided (see flip_geodesic.h the complete listing).

void FlipEdgeNetwork::rewind()

Undo any flips which have been performed, resetting the data structure to its initial state.

Note that you must set FlipEdgeNetwork::supportRewinding = true right after constructure (before doing any shortening) to be able to use rewind().

bool FlipEdgeNetwork::edgeInPath(Edge e)

Test whether a given edge from the intrinsic triangulation appears in the path.

bool FlipEdgeNetwork::halfedgeInPath(Halfedge he)

Test whether a given halfedge from the intrinsic triangulation appears in the path.

double FlipEdgeNetwork::length()

Measure the total length of the path network.

double FlipEdgeNetwork::minAngleIsotopy()

Measure the smallest angle in the path (not blocked by isotopy, aka the constraint that curves do not cross).

void FlipEdgeNetwork::savePathOBJLine(std::string filenamePrefix, bool withAll = false)

Write the path to file as line entries in an obj file.

Citation

If these routines contribute to academic work, please cite:

@article{sharp2020you,
  title={You can find geodesic paths in triangle meshes by just flipping edges},
  author={Sharp, Nicholas and Crane, Keenan},
  journal={ACM Transactions on Graphics (TOG)},
  volume={39},
  number={6},
  pages={1--15},
  year={2020},
  publisher={ACM New York, NY, USA}
}