diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt
index b0958a8ff4f5..9f95dc2cb98a 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt
@@ -152,3 +152,14 @@ if(TARGET CGAL::Ceres_support)
 else()
   message(STATUS "NOTICE: The example 'mesh_smoothing_example' uses the Ceres library, and will not be compiled.")
 endif()
+
+
+find_path(CYCODEBASE_INCLUDE "cySampleElim.h" DOC "Path to the cyCodeBase directory (https://github.com/cemyuksel/cyCodeBase)")
+
+if(NOT CYCODEBASE_INCLUDE)
+  message( STATUS "poisson_eliminate needs the cyCodeBase library (https://github.com/cemyuksel/cyCodeBase)")
+else()
+  create_single_source_cgal_program("poisson_eliminate.cpp")
+  target_compile_definitions( poisson_eliminate CGAL_USE_CY )
+  target_include_directories( poisson_eliminate AFTER PRIVATE ${CYCODEBASE_INCLUDE} )
+endif()
\ No newline at end of file
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/poisson_eliminate.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/poisson_eliminate.cpp
new file mode 100644
index 000000000000..89a8b51359f9
--- /dev/null
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/poisson_eliminate.cpp
@@ -0,0 +1,35 @@
+#include <vector>
+#include <iostream>
+#include <fstream>
+#include <filesystem>
+#include <string>
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/IO/polygon_mesh_io.h>
+#include <CGAL/Polygon_mesh_processing/distance.h>
+#include <CGAL/Polygon_mesh_processing/poisson_eliminate.h>
+
+typedef CGAL::Simple_cartesian<double> K;
+typedef K::Point_3 Point_3;
+typedef CGAL::Surface_mesh<Point_3> Mesh;
+
+
+int main(int argc, char* argv[])
+{
+  std::string filename = std::filesystem::path(argv[1]).stem().string();
+  Mesh sm;
+  CGAL::IO::read_polygon_mesh(argv[1], sm);
+
+  std::vector<Point_3> points;
+
+  CGAL::Polygon_mesh_processing::poisson_eliminate(sm, std::back_inserter(points));
+
+  std::string poisson_points = filename+"-poisson.xyz";
+  std::ofstream out(poisson_points);
+  out.precision(17);
+  for(const Point_3& p : points){
+    out << p << std::endl;
+  }
+  return 0;
+}
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/poisson_eliminate.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/poisson_eliminate.h
new file mode 100644
index 000000000000..f7e7f9595168
--- /dev/null
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/poisson_eliminate.h
@@ -0,0 +1,80 @@
+// Copyright (c) 2024 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)     : Andreas Fabri
+
+#ifndef CGAL_POLYGON_MESH_PROCESSING_POISSON_ELIMINATE_H
+#define CGAL_POLYGON_MESH_PROCESSING_POISSON_ELIMINATE_H
+
+#include <CGAL/license/Polygon_mesh_processing/distance.h>
+
+#ifdef CGAL_USE_CY
+
+#include <cyVector.h>
+#include <cySampleElim.h>
+
+
+#include <CGAL/Polygon_mesh_processing/distance.h>
+#include <CGAL/Polygon_mesh_processing/measure.h>
+
+
+namespace CGAL {
+namespace Polygon_mesh_processing {
+
+template <class TriangleMesh, class OutputIterator, class NamedParameters = parameters::Default_named_parameters>
+poisson_eliminate(const TriangleMesh& sm, OutputIterator out, const NamedParameters& np = parameters::default_values())
+{
+  typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type             GeomTraits;
+  typedef typename GeomTraits::Point_3                                            Point_3;
+
+  Bbox_3 bb = bbox_3(sm.points().begin(), sm.points().end());
+  cy::Vec3d bl(bb.xmin(), bb.ymin(), bb.zmin());
+  cy::Vec3d tr(bb.xmax(), bb.ymax(), bb.zmax());
+
+  std::vector<cy::Vec3d> inputPoints, outputPoints;
+  std::vector<Point_3> points;
+
+  // @todo write with a transform_iterator directly into inputPoints
+  sample_triangle_mesh(sm,
+                       std::back_inserter(points),
+                       CGAL::parameters::number_of_points_on_faces(2* num_vertices(sm))
+                         .do_sample_vertices(false)
+                         .do_sample_edges(false));
+  double area = CGAL::Polygon_mesh_processing::area(sm);
+
+  for(int i = 0; i < points.size(); ++i){
+    inputPoints.push_back(cy::Vec3d(to_double(points[i].x()), to_double(points[i].y()), to_double(points[i].z())));
+  }
+
+  outputPoints.resize(num_vertices(sm)/2);
+
+  cy::WeightedSampleElimination< cy::Vec3d, double, 3, int > wse;
+  wse.SetBoundsMin(bl);
+  wse.SetBoundsMax(tr);
+  bool isProgressive = true;
+
+  double d_max = 2 * wse.GetMaxPoissonDiskRadius( 2, outputPoints.size(), area );
+
+  wse.Eliminate( inputPoints.data(), inputPoints.size(),
+                 outputPoints.data(), outputPoints.size(),
+                 isProgressive,
+                 d_max, 2 );
+
+  for (const cy::Vec3d& p : outputPoints){
+    *out++ = Point_3(p.x, p.y, p.z);
+  }
+}
+
+} // namespace Polygon_mesh_processing
+} // namespace CGAL
+
+#endif // ifdef CGAL_USE_CY
+
+#endif // CGAL_POLYGON_MESH_PROCESSING_POISSON_ELIMINATE_H