BGL - add CGAL::shortest_path(vs, vt, mesh)

BGL/doc/BGL/Doxyfile.in

@@ -19,7 +19,8 @@ INPUT += ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/IO/polygon_mesh_io.h \
          ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/boost/graph/METIS/partition_graph.h \
          ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/boost/graph/METIS/partition_dual_graph.h \
          ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/boost/graph/alpha_expansion_graphcut.h \
-         ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/boost/graph/graph_traits_inheritance_macros.h
+         ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/boost/graph/graph_traits_inheritance_macros.h \
+         ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/boost/graph/shortest_path.h
 
 
 EXAMPLE_PATH               +=  ${CGAL_Surface_mesh_skeletonization_EXAMPLE_DIR} \

BGL/doc/BGL/PackageDescription.txt

@@ -467,6 +467,9 @@ the requirement for traversal of all faces in a graph.
 /// \defgroup PkgBGLPartition Partitioning Operations
 /// \ingroup PkgBGLRef
 
+/// \defgroup PkgBGLTraversal Graph Traversal
+/// \ingroup PkgBGLRef
+
 /// \defgroup PkgBGLIOFct I/O Functions
 /// \ingroup PkgBGLRef
 
@@ -586,6 +589,12 @@ Methods to split a mesh into subdomains, using the library
 implementation.
 */
 
+/*!
+\addtogroup PkgBGLTraversal
+
+Methods to traverse a graph, for example to find the shortest path between two vertices.
+*/
+
 /*!
 \addtogroup PkgBGLIOFct
 
@@ -761,6 +770,9 @@ user might encounter.
 \cgalCRPSection{Conversion Functions}
 - `CGAL::split_graph_into_polylines()`
 
+\cgalCRPSection{Graph Traversal}
+- `CGAL::shortest_path_between_two_vertices()`
+
 \cgalCRPSection{Graph Adaptors}
 - `CGAL::Dual`
 - `CGAL::Face_filtered_graph`

BGL/examples/BGL_surface_mesh/CMakeLists.txt

@@ -9,6 +9,7 @@ create_single_source_cgal_program("seam_mesh.cpp")
 create_single_source_cgal_program("write_inp.cpp")
 create_single_source_cgal_program("surface_mesh_dual.cpp")
 create_single_source_cgal_program("connected_components.cpp")
+create_single_source_cgal_program("shortest_path.cpp")
 
 find_package(METIS QUIET)
 include(CGAL_METIS_support)

BGL/examples/BGL_surface_mesh/shortest_path.cpp

@@ -0,0 +1,80 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/boost/graph/graph_traits_Surface_mesh.h>
+
+#include <CGAL/boost/graph/shortest_path.h>
+#include <CGAL/IO/polygon_mesh_io.h>
+
+
+#include <string>
+#include <vector>
+#include <fstream>
+#include <exception>
+#include <algorithm>
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Point = K::Point_3;
+using Mesh = CGAL::Surface_mesh<Point>;
+
+using vertex_descriptor = boost::graph_traits<Mesh>::vertex_descriptor;
+using edge_descriptor = boost::graph_traits<Mesh>::edge_descriptor;
+using halfedge_descriptor = boost::graph_traits<Mesh>::halfedge_descriptor;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+namespace params = CGAL::parameters;
+
+
+// Example main
+int main(int argc, char** argv)
+{
+  const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/elephant.off");
+
+  // Try building a surface_mesh
+  Mesh sm;
+  bool ok = CGAL::IO::read_polygon_mesh(filename, sm);
+  if (!ok || !sm.is_valid() || sm.is_empty())
+  {
+    std::cerr << "Error: Invalid facegraph" << std::endl;
+    std::cerr << "Filename = " << filename << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  const std::size_t i0 = 0;
+  const std::size_t i1 = num_vertices(sm) / 2;
+
+
+  // Get the vertex descriptors of the source and target vertices
+  const vertex_descriptor vs = *vertices(sm).first;
+  vertex_descriptor vt;
+  std::size_t vid = 0;
+  for (const vertex_descriptor v : vertices(sm))
+  {
+    if (vid++ == i1)
+    {
+      vt = v;
+      break;
+    }
+  }
+
+  std::vector<halfedge_descriptor> halfedge_sequence;
+  CGAL::shortest_path_between_two_vertices(vs, vt, sm,
+    std::back_inserter(halfedge_sequence));
+
+  // dump
+  std::cout << "Shortest path between vertices " << i0 << " and " << i1
+            << " is made of " << halfedge_sequence.size() << " halfedges." << std::endl;
+
+  // Get the property map of the points of the mesh
+  auto vpmap = get(CGAL::vertex_point, sm);
+
+  std::ofstream out("shortest_path.polylines.txt");
+  for (const halfedge_descriptor he : halfedge_sequence)
+  {
+    const vertex_descriptor v0 = source(he, sm);
+    const vertex_descriptor v1 = target(he, sm);
+    out << "2 " << get(vpmap, v0) << " " << get(vpmap, v1) << std::endl;
+  }
+  out.close();
+
+  return EXIT_SUCCESS;
+}

BGL/include/CGAL/boost/graph/shortest_path.h

@@ -0,0 +1,190 @@
+// Copyright (c) 2025  GeometryFactory (France). All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org)
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s)     : Jane Tournois, Andreas Fabri
+//
+
+#ifndef CGAL_BOOST_GRAPH_SHORTEST_PATH_H
+#define CGAL_BOOST_GRAPH_SHORTEST_PATH_H
+
+#include <boost/graph/dijkstra_shortest_paths.hpp>
+
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+
+#include <boost/property_map/property_map.hpp>
+#include <CGAL/boost/graph/properties.h>
+
+#include <vector>
+#include <unordered_map>
+
+namespace CGAL {
+namespace internal {
+
+  /// An exception used while catching a throw that stops Dijkstra's algorithm
+  /// once the shortest path to a target has been found.
+  class Dijkstra_end_exception : public std::exception
+  {
+    const char* what() const throw ()
+    {
+      return "Dijkstra shortest path: reached the target vertex.";
+    }
+  };
+
+  /// Visitor to stop Dijkstra's algorithm once the given target turns 'BLACK',
+  /// that is when the target has been examined through all its incident edges and
+  /// the shortest path is thus known.
+  template<typename Graph, typename VertexEdgeMap>
+  class Stop_at_target_Dijkstra_visitor : boost::default_dijkstra_visitor
+  {
+    using vertex_descriptor = typename boost::graph_traits<Graph>::vertex_descriptor;
+    using edge_descriptor = typename boost::graph_traits<Graph>::edge_descriptor;
+
+  public:
+    vertex_descriptor destination_vd;
+    VertexEdgeMap& relaxed_edges;
+
+    Stop_at_target_Dijkstra_visitor(vertex_descriptor destination_vd,
+                                    VertexEdgeMap& relaxed_edges)
+      : destination_vd(destination_vd), relaxed_edges(relaxed_edges)
+    {}
+
+    void initialize_vertex(const vertex_descriptor& /*s*/, const Graph& /*g*/) const {}
+    void examine_vertex(const vertex_descriptor& /*s*/, const Graph& /*g*/) const {}
+    void examine_edge(const edge_descriptor& /*e*/, const Graph& /*g*/) const {}
+    void edge_relaxed(const edge_descriptor& e, const Graph& g) const
+    {
+      relaxed_edges[target(e, g)] = e;
+    }
+    void discover_vertex(const vertex_descriptor& /*s*/, const Graph& /*g*/) const {}
+    void edge_not_relaxed(const edge_descriptor& /*e*/, const Graph& /*g*/) const {}
+    void finish_vertex(const vertex_descriptor& vd, const Graph& /* g*/) const
+    {
+      if (vd == destination_vd)
+        throw Dijkstra_end_exception();
+    }
+  };
+} // namespace internal
+
+/*!
+* \ingroup PkgBGLTraversal
+* Computes the shortest path between two vertices in a graph `g`.
+* The vertices must belong to the same connected component of `g`.
+*
+* @tparam Graph a model of the concept `HalfedgeListGraph`
+* @tparam OutputIterator an output iterator with value type `boost::graph_traits<Graph>::%halfedge_descriptor`
+* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+*
+* @param vs source vertex
+* @param vt target vertex
+* @param g the graph
+* @param halfedge_sequence_oit the output iterator holding the output sequence
+* of halfedges that form the shortest path from `vs` to `vt` on `g`
+* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+*
+*  \cgalNamedParamsBegin
+*   \cgalParamNBegin{edge_weight_map}
+*    \cgalParamDescription{a property map associating to each edge in the graph its weight or ``length''.
+*                          The weights must all be non-negative.}
+*    \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<PolygonMesh>::%edge_descriptor`
+*                   as key type and `FT` as value type.}
+*    \cgalParamDefault{`get(boost::edge_weight, mesh)`}
+*   \cgalParamNEnd
+*
+*   \cgalParamNBegin{vertex_index_map}
+*     \cgalParamDescription{a property map associating to each vertex of `pmesh` a unique index between `0` and `num_vertices(pmesh) - 1`}
+*     \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<PolygonMesh>::%vertex_descriptor`
+*                    as key type and `std::size_t` as value type}
+*     \cgalParamDefault{an automatically indexed internal map}
+*   \cgalParamNEnd
+*  \cgalNamedParamsEnd
+*/
+template<typename Graph,
+         typename OutputIterator,
+         typename NamedParameters = parameters::Default_named_parameters>
+OutputIterator shortest_path_between_two_vertices(
+  const typename boost::graph_traits<Graph>::vertex_descriptor vs,//source
+  const typename boost::graph_traits<Graph>::vertex_descriptor vt,//target
+  const Graph& g,
+  OutputIterator halfedge_sequence_oit,
+  const NamedParameters& np = parameters::default_values())
+{
+  using vertex_descriptor = typename boost::graph_traits<Graph>::vertex_descriptor;
+  using halfedge_descriptor = typename boost::graph_traits<Graph>::halfedge_descriptor;
+  using edge_descriptor = typename boost::graph_traits<Graph>::edge_descriptor;
+
+  using Pred_umap = std::unordered_map<vertex_descriptor, vertex_descriptor>;
+  using Pred_pmap = boost::associative_property_map<Pred_umap>;
+
+  using parameters::get_parameter;
+  using parameters::choose_parameter;
+
+  const auto w_map = choose_parameter(get_parameter(np, internal_np::edge_weight),
+                                      get(boost::edge_weight, g));
+  const auto vim = get_initialized_vertex_index_map(g, np);
+
+  Pred_umap predecessor;
+  Pred_pmap pred_pmap(predecessor);
+
+  using VEMap = std::unordered_map<vertex_descriptor, edge_descriptor>;
+  VEMap relaxed_edges_map;
+  internal::Stop_at_target_Dijkstra_visitor<Graph, VEMap> vis(vt, relaxed_edges_map);
+  try
+  {
+    boost::dijkstra_shortest_paths(g, vs,
+      boost::predecessor_map(pred_pmap)
+      .visitor(vis)
+      .weight_map(w_map)
+      .vertex_index_map(vim));
+  }
+  catch (const internal::Dijkstra_end_exception& ){}
+
+  // Walk back from target to source and collect vertices along the way
+  struct vertex_on_path
+  {
+    vertex_descriptor vertex;
+    bool is_constrained;
+  };
+
+  std::vector<vertex_descriptor> constrained_vertices = { vs };
+  vertex_descriptor t = vt;
+  std::vector<vertex_on_path> path;
+  do
+  {
+    const bool is_new_vertex = (constrained_vertices.end()
+      == std::find(constrained_vertices.begin(), constrained_vertices.end(), t));
+
+    vertex_on_path vop;
+    vop.vertex = t;
+    vop.is_constrained = !is_new_vertex;
+    path.push_back(vop);
+
+    t = get(pred_pmap, t);
+  }
+  while (t != vs);
+
+  // Add the last vertex
+  vertex_on_path vop;
+  vop.vertex = constrained_vertices.back();
+  vop.is_constrained = true;
+  path.push_back(vop);
+
+  // Display path
+  for (auto path_it = path.begin(); path_it != path.end() - 1; ++path_it)
+  {
+    const auto map_it = vis.relaxed_edges.find(path_it->vertex);
+    if (map_it != vis.relaxed_edges.end())
+      *halfedge_sequence_oit++ = halfedge(map_it->second, g);
+  }
+  return halfedge_sequence_oit;
+}
+
+} // namespace CGAL
+
+
+#endif //CGAL_BOOST_GRAPH_SHORTEST_PATH_H