import poisson disk sampling from Brad's branch:

gsoc2024-triangle_mesh_sampling-bmccoy@8db7ae84e5

Comment out all the geodesic part for now

Author: sloriot

Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt

@@ -21,6 +21,7 @@ create_single_source_cgal_program("polyhedral_envelope.cpp" )
 create_single_source_cgal_program("polyhedral_envelope_of_triangle_soup.cpp" )
 create_single_source_cgal_program("polyhedral_envelope_mesh_containment.cpp" )
 create_single_source_cgal_program("sample_example.cpp" )
+create_single_source_cgal_program("poisson_sampling_sm_example.cpp" )
 create_single_source_cgal_program("self_intersections_example.cpp" )
 create_single_source_cgal_program("volume_connected_components.cpp")
 

Polygon_mesh_processing/examples/Polygon_mesh_processing/poisson_sampling_sm_example.cpp

@@ -0,0 +1,50 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polygon_mesh_processing/distance.h>
+
+#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
+#include <iostream>
+#include <string>
+#include <vector>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel   K;
+typedef K::Point_3                                            Point;
+
+typedef CGAL::Surface_mesh<Point>                             Mesh;
+typedef boost::graph_traits<Mesh>::vertex_descriptor          vertex_descriptor;
+typedef boost::graph_traits<Mesh>::face_descriptor            face_descriptor;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+int main(int argc, char* argv[])
+{
+  const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/eight.off");
+
+  Mesh mesh;
+  if(!PMP::IO::read_polygon_mesh(filename, mesh))
+  {
+    std::cerr << "Invalid input." << std::endl;
+    return 1;
+  }
+
+  const double sampling_radius = (argc > 2) ? std::atof(argv[2]) : 0.009;
+
+  const std::size_t number_of_darts = (argc > 3) ? std::atof(argv[3]) : 20;
+
+  CGAL::Random random_seed = (argc > 4) ? std::atof(argv[4]) : CGAL::get_default_random();
+
+  std::vector<Point> points;
+  PMP::sample_triangle_mesh(mesh,
+                            std::back_inserter(points),
+                            CGAL::parameters::use_poisson_disk_sampling_euclidean(true)
+                                             .sampling_radius(sampling_radius)
+                                             .number_of_darts(number_of_darts)
+                                             .random_seed(random_seed));
+
+  std::ofstream out("sampling.xyz");
+  out << std::setprecision(17);
+  std::copy(points.begin(), points.end(), std::ostream_iterator<Point>(out, "\n"));
+
+  std::cout << points.size() << std::endl;
+  return 0;
+}

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h

@@ -19,6 +19,7 @@
 #include <CGAL/Polygon_mesh_processing/internal/AABB_traversal_traits_with_Hausdorff_distance.h>
 #include <CGAL/Polygon_mesh_processing/measure.h>
 #include <CGAL/Polygon_mesh_processing/bbox.h>
+#include <CGAL/Polygon_mesh_processing/internal/poisson_disk_sampling.h>
 
 #include <CGAL/AABB_tree.h>
 #include <CGAL/AABB_traits_3.h>
@@ -476,6 +477,33 @@ struct Triangle_structure_sampler_for_triangle_mesh
     min_sq_edge_length = (std::numeric_limits<double>::max)();
   }
 
+  void procede()
+  {
+    using parameters::choose_parameter;
+    using parameters::get_parameter;
+    using parameters::is_default_parameter;
+
+    bool use_pds_e = choose_parameter(get_parameter(this->np, internal_np::use_poisson_disk_sampling_euclidean), false);
+//    bool use_pds_g = choose_parameter(get_parameter(this->np, internal_np::use_poisson_disk_sampling_geodesic), false);
+    double sampling_radius = choose_parameter(get_parameter(this->np, internal_np::sampling_radius), 1.);
+    std::size_t number_of_darts = choose_parameter(get_parameter(this->np, internal_np::number_of_darts),30.);
+    CGAL::Random random_seed = choose_parameter(get_parameter(this->np, internal_np::random_seed),CGAL::get_default_random());
+
+
+    if (use_pds_e /* || use_pds_g */)
+    {
+      std::vector<typename GeomTraits::Point_3> points = /* use_pds_e ? */
+                    internal::poisson_disk_sampling<GeomTraits, internal::EUCLIDEAN_DISTANCE>(tm, sampling_radius,number_of_darts,random_seed)
+              /*  : internal::poisson_disk_sampling<GeomTraits, internal::GEODESIC_DISTANCE>(tm, sampling_radius,number_of_darts,random_seed) */;
+      std::copy(points.begin(), points.end(), this->out);
+    }
+    else
+    {
+      static_cast<Base*>(this)->procede();
+    }
+  }
+
+
   std::pair<TriangleIterator, TriangleIterator> get_range()
   {
     return std::make_pair(faces(tm).begin(), faces(tm).end());
@@ -890,6 +918,36 @@ struct Triangle_structure_sampler_for_triangle_soup
  *     \cgalParamType{unsigned int}
  *     \cgalParamDefault{`0`}
  *   \cgalParamNEnd
+ *
+ *   \cgalParamNBegin{use_poisson_disk_sampling_euclidean}
+ *     \cgalParamDescription{if `true` is passed, the Euclidean distance is used to compute the distance between sampled points.}
+ *     \cgalParamType{Boolean}
+ *     \cgalParamDefault{`false`}
+ *   \cgalParamNEnd
+ *
+ *   \cgalParamNBegin{use_poisson_disk_sampling_geodesic}
+ *     \cgalParamDescription{if `true` is passed, the approximate geodesic distance is used to compute the distance between
+ *                          sampled points.}
+ *     \cgalParamType{Boolean}
+ *     \cgalParamDefault{`false`}
+ *     \cgalParamExtra{The geodesic distance is approximated using the 'locally_shortest_path'
+ *                     function.}
+ *   \cgalParamNEnd
+ *
+ *   \cgalParamNBegin{sampling_radius}
+ *     \cgalParamDescription{a value used by Poisson disk sampling to specify the minimum allowable distance between
+ *                          points in the sample.}
+ *     \cgalParamType{double}
+ *     \cgalParamDefault{`1`}
+ *   \cgalParamNEnd
+ *
+ *   \cgalParamNBegin{number_of_darts}
+ *     \cgalParamDescription{a value used by Poisson disk sampling to specify the number of attempts to find a point in the annulus
+ *                          around a sample point that is sufficiently far from all other points in the sample.}
+ *     \cgalParamType{std::size_t}
+ *     \cgalParamDefault{`30`}
+ *   \cgalParamNEnd
+ *
  * \cgalNamedParamsEnd
  *
  * @see `CGAL::Polygon_mesh_processing::sample_triangle_soup()`
@@ -970,7 +1028,7 @@ sample_triangle_mesh(const TriangleMesh& tm,
  *                     of the smallest non-null edge of the soup or the value passed to the named parameter
  *                     `grid_spacing`.}
  *   \cgalParamNEnd
- *     \cgalParamNBegin{use_monte_carlo_sampling}
+ *   \cgalParamNBegin{use_monte_carlo_sampling}
  *     \cgalParamDescription{if `true` is passed, points are generated randomly in each triangle.}
  *     \cgalParamType{Boolean}
  *     \cgalParamDefault{`false`}

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/poisson_disk_sampling.h

@@ -0,0 +1,250 @@
+// Copyright (c) 2025 GeometryFactory (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Sebastien Loriot
+//                 Bradley McCoy
+
+#ifndef CGAL_POLYGON_MESH_PROCESSING_INTERNAL_POISSON_DISK_SAMPLING_H
+#define CGAL_POLYGON_MESH_PROCESSING_INTERNAL_POISSON_DISK_SAMPLING_H
+
+
+#include <CGAL/point_generators_3.h>
+#include <CGAL/Polygon_mesh_processing/locate.h>
+#include <CGAL/boost/graph/Face_filtered_graph.h>
+#include <CGAL/squared_distance_3.h>
+#include <CGAL/boost/graph/helpers.h>
+#include <CGAL/Dynamic_property_map.h>
+#ifdef CGAL_USE_BSURF
+#include <CGAL/Polygon_mesh_processing/Bsurf/locally_shortest_path.h>
+#endif
+
+#include <vector>
+#include <queue>
+
+namespace CGAL {
+namespace Polygon_mesh_processing {
+namespace internal {
+
+enum Distance_version { EUCLIDEAN_DISTANCE, GEODESIC_DISTANCE };
+
+template <class GeomTraits>
+double euclideanDistancePoints(const typename GeomTraits::Point_3& source,
+                               const typename GeomTraits::Point_3& target)
+{
+  return sqrt(CGAL::squared_distance(source,target));
+}
+
+#ifdef CGAL_USE_BSURF
+template <class GeomTraits, class TriangleMesh>
+double geodesiceApproximation(const Face_location<TriangleMesh, typename GeomTraits::FT>& source,
+                              const Face_location<TriangleMesh, typename GeomTraits::FT>& target,
+                              const TriangleMesh& mesh,
+                              const Dual_geodesic_solver<double>& solver)
+{
+  std::vector<Edge_location<TriangleMesh, typename GeomTraits::FT>> edge_locations;
+  CGAL::Polygon_mesh_processing::locally_shortest_path<double>(source, target, mesh, edge_locations, solver);
+
+  return path_length<GeomTraits>(edge_locations,source,target,mesh);
+}
+#endif
+
+//function to switch between geodesic and Euclidean distance
+template <class GeomTraits, Distance_version V, class TriangleMesh>
+double distancePoints(const TriangleMesh& /* mesh */,
+                      const typename GeomTraits::Point_3& source,
+                      const typename GeomTraits::Point_3& target,
+                      const Face_location<TriangleMesh, typename GeomTraits::FT>& /* start */,
+                      const Face_location<TriangleMesh, typename GeomTraits::FT>& /* end */
+#ifdef CGAL_USE_BSURF
+                      , const Dual_geodesic_solver<double>& solver
+#endif
+                      )
+{
+#ifdef CGAL_USE_BSURF
+  if constexpr (V==GEODESIC_DISTANCE)
+    return geodesiceApproximation<GeomTraits>(start, end, mesh, solver);
+#endif
+  if constexpr (V==EUCLIDEAN_DISTANCE)
+    return euclideanDistancePoints<GeomTraits>(source, target);
+  return 0;
+}
+
+
+template <class GeomTraits, Distance_version V, class TriangleMesh>
+std::vector<typename boost::graph_traits<TriangleMesh>::face_descriptor>
+faces_in_sub_mesh(const typename GeomTraits::Point_3& c,
+                  typename boost::graph_traits<TriangleMesh>::face_descriptor fc,
+                  const TriangleMesh& mesh,
+                  double minDistance)
+{
+//  using Point = typename GeomTraits::Point_3;
+  using Graph_traits = boost::graph_traits<TriangleMesh>;
+  using face_descriptor = typename Graph_traits::face_descriptor;
+  using halfedge_descriptor = typename Graph_traits::halfedge_descriptor;
+
+  //Begin flooding
+  face_descriptor fd = fc;
+
+  std::vector<bool> selected(num_faces(mesh), false);
+  std::vector<face_descriptor> selection;
+  selected[fd] = true;
+  selection.push_back(fd);
+
+  auto do_queue_edge = [&](halfedge_descriptor h)
+  {
+    halfedge_descriptor hopp=opposite(h, mesh);
+    if (is_border(hopp, mesh) || selected[face(hopp, mesh)]) return false;
+
+    typename GeomTraits::Segment_3 edge(mesh.point(source(h,mesh)), mesh.point(target(h,mesh)));
+
+//      return (distancePoints<V>(mesh, mesh.point(source(h,mesh)), c, tree)< 3*minDistance ||
+//              distancePoints<V>(mesh, mesh.point(target(h,mesh)), c, tree)< 3*minDistance);
+
+
+    return typename GeomTraits::Compare_squared_distance_3()(c, edge, sqrt(3*minDistance))!=CGAL::LARGER;
+  };
+
+
+  std::vector<halfedge_descriptor> queue;
+  for (halfedge_descriptor h : CGAL::halfedges_around_face(halfedge(fd, mesh), mesh))
+    if (do_queue_edge(h))
+      queue.push_back(opposite(h, mesh));
+
+  while (!queue.empty())
+  {
+    halfedge_descriptor h = queue.back();
+    face_descriptor f = face(h, mesh);
+    queue.pop_back();
+    if (!selected[f])
+    {
+      selected[f]=true;
+      selection.push_back(f);
+    }
+
+    h=next(h, mesh);
+    if (do_queue_edge(h)) queue.push_back(opposite(h, mesh));
+    h=next(h, mesh);
+    if (do_queue_edge(h)) queue.push_back(opposite(h, mesh));
+  }
+
+  return selection;
+}
+
+template <class GeomTraits, Distance_version V, class TriangleMesh, class PointsPerFaceMap>
+bool
+is_far_enough(const typename GeomTraits::Point_3 c,
+              const Face_location<TriangleMesh, typename GeomTraits::FT> c_location,
+              const TriangleMesh& mesh,
+              double minDistance,
+              std::vector<typename boost::graph_traits<TriangleMesh>::face_descriptor> selection,
+              const PointsPerFaceMap& face_points
+#ifdef CGAL_USE_BSURF
+              , const Dual_geodesic_solver<double>& solver
+#endif
+              )
+{
+  for (typename boost::graph_traits<TriangleMesh>::face_descriptor f : selection)
+  {
+    for(const typename GeomTraits::Point_3& p  : face_points[f])
+    {
+      Face_location<TriangleMesh, typename GeomTraits::FT> p_location = locate_in_face(p, f, mesh);
+      //Todo: Ask why  is distancePoints for Euclidean so much slower then just
+      //calling euclideanDistancePoints?
+      if (distancePoints<GeomTraits,V>(mesh, c, p, c_location, p_location
+#ifdef CGAL_USE_BSURF
+          , solver
+#endif
+                                      ) < minDistance)
+      //if (euclideanDistancePoints(c, p) < minDistance)
+      //if (geodesiceApproximation(c_location, p_location, mesh) < minDistance)
+          return false;
+    }
+  }
+  return true;
+}
+
+template <class GeomTraits, Distance_version V, class TriangleMesh>
+std::vector<typename GeomTraits::Point_3>
+poisson_disk_sampling(const TriangleMesh& mesh,
+                      double minDistance,
+                      std::size_t kMaxTries,
+                      CGAL::Random& r)
+{
+  using Point = typename GeomTraits::Point_3;
+  using Graph_traits = boost::graph_traits<TriangleMesh>;
+  using face_descriptor = typename Graph_traits::face_descriptor;
+  using FT = typename GeomTraits::FT;
+  using Face_location = Face_location<TriangleMesh, FT>;
+
+#ifdef CGAL_USE_BSURF
+  Dual_geodesic_solver<double> solver;
+  if constexpr (V==GEODESIC_DISTANCE)
+    init_geodesic_dual_solver(solver, mesh);
+#endif
+
+  std::vector<Point> points;
+  std::queue<std::pair<Point, face_descriptor>> activePoints;
+
+  Random_points_in_triangle_mesh_3<TriangleMesh> g(mesh, r);
+  Point c = *g;
+  face_descriptor fc = g.last_item_picked();
+
+  std::vector<std::vector<Point>> face_points(num_faces(mesh));
+
+  face_points[fc].push_back(c);
+  activePoints.push(std::make_pair(c, fc));
+  points.push_back(c);
+
+  Point currentPoint;
+  face_descriptor currentFace;
+  while (!activePoints.empty())
+  {
+    std::tie(currentPoint, currentFace) = activePoints.front();
+    activePoints.pop();
+
+    std::vector<face_descriptor> selection = faces_in_sub_mesh<GeomTraits,V>(currentPoint, currentFace, mesh, minDistance);
+    Face_location current_location = locate_in_face(currentPoint, currentFace, mesh);
+    if (current_location.second[0]<0) current_location.second[0]=0;
+    if (current_location.second[1]<0) current_location.second[1]=0;
+    if (current_location.second[2]<0) current_location.second[2]=0;
+    std::size_t k = 0;
+    while (k < kMaxTries)
+    {
+      double angle = 2 * CGAL_PI * r.get_double();
+      double distance = minDistance + minDistance * r.get_double();
+      typename GeomTraits::Vector_2 dir(cos(angle),sin(angle));
+      #warning find alternative
+      CGAL_USE(distance);
+      std::vector<Face_location> path /* = straightest_geodesic<GeomTraits>(current_location, dir, distance, mesh) */ ;
+      Face_location newLocation = path.back();
+      Point newPoint = construct_point(path.back(), mesh);
+      if(is_far_enough<GeomTraits,V>(newPoint, newLocation, mesh, minDistance, selection, make_random_access_property_map(face_points)
+#ifdef CGAL_USE_BSURF
+        , solver
+#endif
+      ))
+      {
+        face_points[path.back().first].push_back(newPoint);
+        activePoints.push(std::make_pair(newPoint,path.back().first));
+        points.push_back(newPoint);
+      }else
+      {
+        ++k;
+      }
+    }
+  }
+
+  return points;
+}
+
+} } } // CGAL::Polygon_mesh_processing::internal
+
+
+#endif //CGAL_POLYGON_MESH_PROCESSING_INTERNAL_POISSON_DISK_SAMPLING_H