feat(StraightSkeleton): Add partitioning algorithm using Straight Ske…

…leton

src/algorithm/straightSkeleton.cpp

@@ -22,6 +22,8 @@
 #include "SFCGAL/algorithm/tesselate.h"
 #include "SFCGAL/algorithm/translate.h"
 
+#include <CGAL/Arr_segment_traits_2.h>
+#include <CGAL/Arrangement_2.h>
 #include <CGAL/Straight_skeleton_converter_2.h>
 #include <CGAL/Surface_mesh.h>
 #include <CGAL/create_straight_skeleton_from_polygon_with_holes_2.h>
@@ -37,6 +39,8 @@ using Polygon_2            = CGAL::Polygon_2<Kernel>;
 using Polygon_with_holes_2 = CGAL::Polygon_with_holes_2<Kernel>;
 using Straight_skeleton_2  = CGAL::Straight_skeleton_2<Kernel>;
 using Mesh                 = CGAL::Surface_mesh<Point_3>;
+using Segment_2            = Kernel::Segment_2;
+using Arrangement_2 = CGAL::Arrangement_2<CGAL::Arr_segment_traits_2<Kernel>>;
 
 #if CGAL_VERSION_MAJOR < 6
 template <class T>
@@ -443,5 +447,106 @@ extrudeStraightSkeleton(const Geometry &g, double building_height,
       new PolyhedralSurface(building->as<Solid>().exteriorShell())};
   result->addPolygons(*roof);
   return result;
-};
+}
+
+auto
+straightSkeletonPartition(const Geometry &g, bool autoOrientation)
+    -> std::unique_ptr<MultiPolygon>
+{
+  SFCGAL_ASSERT_GEOMETRY_VALIDITY_2D(g);
+
+  std::unique_ptr<MultiPolygon> result(new MultiPolygon);
+
+  switch (g.geometryTypeId()) {
+  case TYPE_TRIANGLE:
+    return straightSkeletonPartition(g.as<Triangle>().toPolygon(),
+                                     autoOrientation);
+  case TYPE_POLYGON:
+    return straightSkeletonPartition(g.as<Polygon>(), autoOrientation);
+  case TYPE_MULTIPOLYGON:
+    return straightSkeletonPartition(g.as<MultiPolygon>(), autoOrientation);
+  default:
+    BOOST_THROW_EXCEPTION(
+        Exception("Geometry must be a Polygon or MultiPolygon"));
+  }
+
+  return result;
+}
+
+auto
+straightSkeletonPartition(const MultiPolygon &g, bool autoOrientation)
+    -> std::unique_ptr<MultiPolygon>
+{
+  std::unique_ptr<MultiPolygon> result(new MultiPolygon);
+
+  for (size_t i = 0; i < g.numGeometries(); i++) {
+    std::unique_ptr<MultiPolygon> partitioned =
+        straightSkeletonPartition(g.polygonN(i), autoOrientation);
+    for (size_t j = 0; j < partitioned->numGeometries(); j++) {
+      result->addGeometry(partitioned->geometryN(j));
+    }
+  }
+
+  return result;
+}
+
+auto
+straightSkeletonPartition(const Polygon &g, bool /*autoOrientation*/)
+    -> std::unique_ptr<MultiPolygon>
+{
+  std::unique_ptr<MultiPolygon> result(new MultiPolygon);
+
+  if (g.isEmpty()) {
+    return result;
+  }
+
+  Kernel::Vector_2                      trans;
+  Polygon_with_holes_2 const            polygon  = preparePolygon(g, trans);
+  SHARED_PTR<Straight_skeleton_2> const skeleton = straightSkeleton(polygon);
+
+  if (skeleton == nullptr) {
+    BOOST_THROW_EXCEPTION(Exception("CGAL failed to create straightSkeleton"));
+  }
+
+  // Function to create a polygon from a face
+  auto create_polygon_from_face =
+      [&trans](const Straight_skeleton_2::Face_const_handle &face) {
+        std::vector<Point>                         points;
+        Straight_skeleton_2::Halfedge_const_handle start   = face->halfedge();
+        Straight_skeleton_2::Halfedge_const_handle current = start;
+        do {
+          const Point_2 &p = current->vertex()->point();
+          points.emplace_back(CGAL::to_double(p.x()), CGAL::to_double(p.y()));
+          current = current->next();
+        } while (current != start);
+
+        SFCGAL::Polygon poly(SFCGAL::LineString(std::move(points)));
+        algorithm::translate(poly, trans);
+        return poly;
+      };
+
+  // Function to check if a face corresponds to a hole
+  auto is_hole_face = [&polygon](
+                          const Straight_skeleton_2::Face_const_handle &face) {
+    Point_2 test_point = face->halfedge()->vertex()->point();
+    for (auto hit = polygon.holes_begin(); hit != polygon.holes_end(); ++hit) {
+      if (hit->bounded_side(test_point) == CGAL::ON_BOUNDED_SIDE) {
+        return true;
+      }
+    }
+    return false;
+  };
+
+  // Iterate through the faces of the skeleton
+  for (auto face = skeleton->faces_begin(); face != skeleton->faces_end();
+       ++face) {
+    // Skip the faces that correspond to holes
+    if (is_hole_face(face))
+      continue;
+
+    result->addGeometry(create_polygon_from_face(face));
+  }
+
+  return result;
+}
 } // namespace SFCGAL::algorithm

src/algorithm/straightSkeleton.h

@@ -104,6 +104,53 @@ extrudeStraightSkeleton(const Geometry &g, double building_height,
                         double roof_height)
     -> std::unique_ptr<PolyhedralSurface>;
 
+/**
+ * @brief Build a 2D straight skeleton partition for a Geometry
+ * @ingroup public_api
+ * @param[in] g The input geometry
+ * @param[in] autoOrientation Check and fix polygon orientation
+ * @return A unique pointer to a MultiPolygon representing the partitioned
+ * geometry
+ * @throws Exception If CGAL fails to create the straight skeleton
+ * @note Only Triangle, Polygon, and MultiPolygon geometries are supported
+ *
+ * This function creates a partition of the input geometry based on its straight
+ * skeleton. For unsupported geometry types, an empty MultiPolygon is returned.
+ */
+SFCGAL_API std::unique_ptr<MultiPolygon>
+straightSkeletonPartition(const Geometry &g, bool autoOrientation = true);
+
+/**
+ * @brief Build a 2D straight skeleton partition for a Polygon
+ * @ingroup detail
+ * @param[in] g The input polygon
+ * @param[in] autoOrientation Check and fix polygon orientation (not used in
+ * this implementation)
+ * @return A unique pointer to a MultiPolygon representing the partitioned
+ * polygon
+ * @throws Exception If CGAL fails to create the straight skeleton
+ *
+ * This function creates a partition of the input polygon based on its straight
+ * skeleton. It uses CGAL's Arrangement_2 to handle the intersection of skeleton
+ * and polygon edges.
+ */
+SFCGAL_API std::unique_ptr<MultiPolygon>
+straightSkeletonPartition(const Polygon &g, bool autoOrientation = true);
+
+/**
+ * @brief Build a 2D straight skeleton partition for a MultiPolygon
+ * @ingroup detail
+ * @param[in] g The input multi-polygon
+ * @param[in] autoOrientation Check and fix polygon orientation
+ * @return A unique pointer to a MultiPolygon representing the partitioned
+ * multi-polygon
+ *
+ * This function applies the straight skeleton partition to each polygon in the
+ * input multi-polygon and combines the results into a single MultiPolygon.
+ */
+SFCGAL_API std::unique_ptr<MultiPolygon>
+straightSkeletonPartition(const MultiPolygon &g, bool autoOrientation = true);
+
 } // namespace algorithm
 } // namespace SFCGAL
 

src/capi/sfcgal_c.cpp

@@ -1716,3 +1716,23 @@ sfcgal_geometry_rotate_z(const sfcgal_geometry_t *geom, double angle)
       std::unique_ptr<SFCGAL::Geometry> result(g.clone());
       SFCGAL::algorithm::rotateZ(*result, angle); return result.release();)
 }
+
+extern "C" auto
+sfcgal_geometry_straight_skeleton_partition(const sfcgal_geometry_t *geom,
+                                            bool autoOrientation)
+    -> sfcgal_geometry_t *
+{
+  std::unique_ptr<SFCGAL::Geometry> result;
+  try {
+    result = SFCGAL::algorithm::straightSkeletonPartition(
+        *(const SFCGAL::Geometry *)(geom));
+  } catch (std::exception &e) {
+    SFCGAL_WARNING("During straight_skeleton_partition (A, %g) :",
+                   autoOrientation);
+    SFCGAL_WARNING("  with A: %s",
+                   ((const SFCGAL::Geometry *)(geom))->asText().c_str());
+    SFCGAL_ERROR("%s", e.what());
+    return nullptr;
+  }
+  return result.release();
+}

src/capi/sfcgal_c.h

@@ -1065,6 +1065,16 @@ sfcgal_geometry_extrude_polygon_straight_skeleton(const sfcgal_geometry_t *geom,
 SFCGAL_API sfcgal_geometry_t *
 sfcgal_geometry_approximate_medial_axis(const sfcgal_geometry_t *geom);
 
+/**
+ * Returns the straight skeleton partition for the given Polygon
+ * @pre isValid(geom) == true
+ * @post isValid(return) == true
+ * @ingroup capi
+ */
+SFCGAL_API sfcgal_geometry_t *
+sfcgal_geometry_straight_skeleton_partition(const sfcgal_geometry_t *geom,
+                                            bool autoOrientation);
+
 /**
  * Tests the coverage of geom1 and geom2
  * @pre isValid(geom1) == true

test/unit/SFCGAL/algorithm/StraightSkeletonTest.cpp

@@ -417,4 +417,97 @@ BOOST_AUTO_TEST_CASE(testExtrudeStraightSkeletonGenerateBuilding)
   BOOST_CHECK_EQUAL(out->asText(2), expectedWKT);
 }
 
+BOOST_AUTO_TEST_CASE(testStraightSkeletonPartitionLShapedPolygon)
+{
+
+  std::unique_ptr<Geometry> const g(
+      io::readWkt("POLYGON((0 0, 0 2, 1 2, 1 1, 2 1, 2 0, 0 0))"));
+  std::unique_ptr<Geometry> out(algorithm::straightSkeletonPartition(*g));
+  std::string const         expectedWKT(
+      "MULTIPOLYGON(((0.00 0.00,0.50 0.50,0.50 1.50,0.00 2.00)),((2.00 "
+              "0.00,1.50 0.50,0.50 0.50,0.00 0.00)),((2.00 1.00,1.50 0.50,2.00 "
+              "0.00)),((1.00 1.00,0.50 0.50,1.50 0.50,2.00 1.00)),((1.00 2.00,0.50 "
+              "1.50,0.50 0.50,1.00 1.00)),((0.00 2.00,0.50 1.50,1.00 2.00)))");
+  BOOST_CHECK_EQUAL(out->asText(2), expectedWKT);
+}
+
+BOOST_AUTO_TEST_CASE(testStraightSkeletonPartitionSimpleRectangle)
+{
+  std::unique_ptr<Geometry> g(
+      io::readWkt("POLYGON((0 0, 0 2, 3 2, 3 0, 0 0))"));
+  std::unique_ptr<Geometry> out(algorithm::straightSkeletonPartition(*g));
+  BOOST_CHECK(out->is<MultiPolygon>());
+  BOOST_CHECK_EQUAL(out->as<MultiPolygon>().numGeometries(), 4);
+  std::string const expectedWKT(
+      "MULTIPOLYGON(((0.00 0.00,1.00 1.00,0.00 2.00)),((3.00 0.00,2.00 "
+      "1.00,1.00 1.00,0.00 0.00)),((3.00 2.00,2.00 1.00,3.00 0.00)),((0.00 "
+      "2.00,1.00 1.00,2.00 1.00,3.00 2.00)))");
+  BOOST_CHECK_EQUAL(out->asText(2), expectedWKT);
+}
+
+BOOST_AUTO_TEST_CASE(testStraightSkeletonPartitionComplexPolygon)
+{
+  std::unique_ptr<Geometry> g(
+      io::readWkt("POLYGON((0 0, 0 3, 1 3, 1 2, 2 2, 2 3, 3 3, 3 0, 0 0))"));
+  std::unique_ptr<Geometry> out(algorithm::straightSkeletonPartition(*g));
+  BOOST_CHECK(out->is<MultiPolygon>());
+  BOOST_CHECK_EQUAL(out->as<MultiPolygon>().numGeometries(), 8);
+  std::string const expectedWKT(
+      "MULTIPOLYGON(((0.00 0.00,1.00 1.00,0.50 1.50,0.50 2.50,0.00 "
+      "3.00)),((3.00 0.00,2.00 1.00,1.00 1.00,0.00 0.00)),((3.00 3.00,2.50 "
+      "2.50,2.50 1.50,2.00 1.00,3.00 0.00)),((2.00 3.00,2.50 2.50,3.00 "
+      "3.00)),((2.00 2.00,2.50 1.50,2.50 2.50,2.00 3.00)),((1.00 2.00,0.50 "
+      "1.50,1.00 1.00,2.00 1.00,2.50 1.50,2.00 2.00)),((1.00 3.00,0.50 "
+      "2.50,0.50 1.50,1.00 2.00)),((0.00 3.00,0.50 2.50,1.00 3.00)))");
+  BOOST_CHECK_EQUAL(out->asText(2), expectedWKT);
+}
+
+BOOST_AUTO_TEST_CASE(testStraightSkeletonPartitionPolygonWithHole)
+{
+  std::unique_ptr<Geometry> g(io::readWkt(
+      "POLYGON((0 0, 0 4, 4 4, 4 0, 0 0), (1 1, 3 1, 3 3, 1 3, 1 1))"));
+  std::unique_ptr<Geometry> out(algorithm::straightSkeletonPartition(*g));
+  BOOST_CHECK(out->is<MultiPolygon>());
+  std::string const expectedWKT(
+      "MULTIPOLYGON(((0.00 0.00,0.50 0.50,0.50 3.50,0.00 4.00)),((4.00 "
+      "0.00,3.50 0.50,0.50 0.50,0.00 0.00)),((4.00 4.00,3.50 3.50,3.50 "
+      "0.50,4.00 0.00)),((0.00 4.00,0.50 3.50,3.50 3.50,4.00 4.00)),((1.00 "
+      "1.00,0.50 0.50,3.50 0.50,3.00 1.00)),((1.00 3.00,0.50 3.50,0.50 "
+      "0.50,1.00 1.00)),((3.00 3.00,3.50 3.50,0.50 3.50,1.00 3.00)),((3.00 "
+      "1.00,3.50 0.50,3.50 3.50,3.00 3.00)))");
+  BOOST_CHECK_EQUAL(out->asText(2), expectedWKT);
+}
+
+BOOST_AUTO_TEST_CASE(testStraightSkeletonPartitionMultiPolygon)
+{
+  std::unique_ptr<Geometry> g(
+      io::readWkt("MULTIPOLYGON(((0 0, 0 2, 2 2, 2 0, 0 0)), ((3 3, 3 5, 5 5, "
+                  "5 3, 3 3)))"));
+  std::unique_ptr<Geometry> out(algorithm::straightSkeletonPartition(*g));
+  BOOST_CHECK(out->is<MultiPolygon>());
+  std::string const expectedWKT(
+      "MULTIPOLYGON(((0.00 0.00,1.00 1.00,0.00 2.00)),((2.00 0.00,1.00 "
+      "1.00,0.00 0.00)),((2.00 2.00,1.00 1.00,2.00 0.00)),((0.00 2.00,1.00 "
+      "1.00,2.00 2.00)),((3.00 3.00,4.00 4.00,3.00 5.00)),((5.00 3.00,4.00 "
+      "4.00,3.00 3.00)),((5.00 5.00,4.00 4.00,5.00 3.00)),((3.00 5.00,4.00 "
+      "4.00,5.00 5.00)))");
+  BOOST_CHECK_EQUAL(out->asText(2), expectedWKT);
+}
+
+BOOST_AUTO_TEST_CASE(testStraightSkeletonPartitionEmptyPolygon)
+{
+  std::unique_ptr<Geometry> g(io::readWkt("POLYGON EMPTY"));
+  std::unique_ptr<Geometry> out(algorithm::straightSkeletonPartition(*g));
+  BOOST_CHECK(out->is<MultiPolygon>());
+  BOOST_CHECK(out->isEmpty());
+  std::string const expectedWKT("MULTIPOLYGON EMPTY");
+  BOOST_CHECK_EQUAL(out->asText(2), expectedWKT);
+}
+
+BOOST_AUTO_TEST_CASE(testStraightSkeletonPartitionNonPolygonGeometry)
+{
+  std::unique_ptr<Geometry> g(io::readWkt("LINESTRING(0 0, 1 1, 2 2)"));
+  BOOST_CHECK_THROW(algorithm::straightSkeletonPartition(*g), std::exception);
+}
+
 BOOST_AUTO_TEST_SUITE_END()

test/unit/SFCGAL/capi/sfcgal_cTest.cpp

@@ -423,4 +423,28 @@ BOOST_AUTO_TEST_CASE(testRotateX)
   sfcgal_geometry_delete(rotated);
 }
 
+BOOST_AUTO_TEST_CASE(testStraightSkeletonPartitionC)
+{
+
+  sfcgal_set_error_handlers(printf, on_error);
+  std::unique_ptr<Geometry> const geom(
+      io::readWkt("POLYGON((0 0, 0 2, 1 2, 1 1, 2 1, 2 0, 0 0))"));
+  std::string const expectedWKT(
+      "MULTIPOLYGON(((0.00 0.00,0.50 0.50,0.50 1.50,0.00 2.00)),((2.00 "
+      "0.00,1.50 0.50,0.50 0.50,0.00 0.00)),((2.00 1.00,1.50 0.50,2.00 "
+      "0.00)),((1.00 1.00,0.50 0.50,1.50 0.50,2.00 1.00)),((1.00 2.00,0.50 "
+      "1.50,0.50 0.50,1.00 1.00)),((0.00 2.00,0.50 1.50,1.00 2.00)))");
+  sfcgal_geometry_t *result =
+      sfcgal_geometry_straight_skeleton_partition(geom.get(), true);
+  // retrieve wkb from C api
+  char  *wkbApi;
+  size_t wkbLen;
+  sfcgal_geometry_as_text_decim(result, 2, &wkbApi, &wkbLen);
+  std::string strApi(wkbApi, wkbLen);
+
+  // check
+  BOOST_CHECK_EQUAL(expectedWKT, strApi);
+  delete[] wkbApi;
+}
+
 BOOST_AUTO_TEST_SUITE_END()