From d5d2ea3e6c6cec2ae8e77ccd046a63dcc0845cf8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Bartoletti?= Date: Tue, 10 Sep 2024 14:55:31 +0200 Subject: [PATCH] feat(StraightSkeleton): Add partitioning algorithm using Straight Skeleton --- src/algorithm/straightSkeleton.cpp | 107 +++++++++++++++++- src/algorithm/straightSkeleton.h | 47 ++++++++ src/capi/sfcgal_c.cpp | 20 ++++ src/capi/sfcgal_c.h | 10 ++ .../SFCGAL/algorithm/StraightSkeletonTest.cpp | 93 +++++++++++++++ test/unit/SFCGAL/capi/sfcgal_cTest.cpp | 25 ++++ 6 files changed, 301 insertions(+), 1 deletion(-) diff --git a/src/algorithm/straightSkeleton.cpp b/src/algorithm/straightSkeleton.cpp index 6e47400c..beb1e038 100644 --- a/src/algorithm/straightSkeleton.cpp +++ b/src/algorithm/straightSkeleton.cpp @@ -22,6 +22,8 @@ #include "SFCGAL/algorithm/tesselate.h" #include "SFCGAL/algorithm/translate.h" +#include +#include #include #include #include @@ -37,6 +39,8 @@ using Polygon_2 = CGAL::Polygon_2; using Polygon_with_holes_2 = CGAL::Polygon_with_holes_2; using Straight_skeleton_2 = CGAL::Straight_skeleton_2; using Mesh = CGAL::Surface_mesh; +using Segment_2 = Kernel::Segment_2; +using Arrangement_2 = CGAL::Arrangement_2>; #if CGAL_VERSION_MAJOR < 6 template @@ -443,5 +447,106 @@ extrudeStraightSkeleton(const Geometry &g, double building_height, new PolyhedralSurface(building->as().exteriorShell())}; result->addPolygons(*roof); return result; -}; +} + +auto +straightSkeletonPartition(const Geometry &g, bool autoOrientation) + -> std::unique_ptr +{ + SFCGAL_ASSERT_GEOMETRY_VALIDITY_2D(g); + + std::unique_ptr result(new MultiPolygon); + + switch (g.geometryTypeId()) { + case TYPE_TRIANGLE: + return straightSkeletonPartition(g.as().toPolygon(), + autoOrientation); + case TYPE_POLYGON: + return straightSkeletonPartition(g.as(), autoOrientation); + case TYPE_MULTIPOLYGON: + return straightSkeletonPartition(g.as(), 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 +{ + std::unique_ptr result(new MultiPolygon); + + for (size_t i = 0; i < g.numGeometries(); i++) { + std::unique_ptr 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 +{ + std::unique_ptr result(new MultiPolygon); + + if (g.isEmpty()) { + return result; + } + + Kernel::Vector_2 trans; + Polygon_with_holes_2 const polygon = preparePolygon(g, trans); + SHARED_PTR 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 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 diff --git a/src/algorithm/straightSkeleton.h b/src/algorithm/straightSkeleton.h index ec170df7..4012c0da 100644 --- a/src/algorithm/straightSkeleton.h +++ b/src/algorithm/straightSkeleton.h @@ -104,6 +104,53 @@ extrudeStraightSkeleton(const Geometry &g, double building_height, double roof_height) -> std::unique_ptr; +/** + * @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 +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 +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 +straightSkeletonPartition(const MultiPolygon &g, bool autoOrientation = true); + } // namespace algorithm } // namespace SFCGAL diff --git a/src/capi/sfcgal_c.cpp b/src/capi/sfcgal_c.cpp index 63effbd5..7df8a8c2 100644 --- a/src/capi/sfcgal_c.cpp +++ b/src/capi/sfcgal_c.cpp @@ -1535,3 +1535,23 @@ sfcgal_geometry_visibility_segment(const sfcgal_geometry_t *polygon, return result.release(); } + +extern "C" auto +sfcgal_geometry_straight_skeleton_partition(const sfcgal_geometry_t *geom, + bool autoOrientation) + -> sfcgal_geometry_t * +{ + std::unique_ptr 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(); +} diff --git a/src/capi/sfcgal_c.h b/src/capi/sfcgal_c.h index 7c9b71ca..1abb4871 100644 --- a/src/capi/sfcgal_c.h +++ b/src/capi/sfcgal_c.h @@ -1058,6 +1058,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 diff --git a/test/unit/SFCGAL/algorithm/StraightSkeletonTest.cpp b/test/unit/SFCGAL/algorithm/StraightSkeletonTest.cpp index 3956a5ce..1558a65b 100644 --- a/test/unit/SFCGAL/algorithm/StraightSkeletonTest.cpp +++ b/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 const g( + io::readWkt("POLYGON((0 0, 0 2, 1 2, 1 1, 2 1, 2 0, 0 0))")); + std::unique_ptr 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 g( + io::readWkt("POLYGON((0 0, 0 2, 3 2, 3 0, 0 0))")); + std::unique_ptr out(algorithm::straightSkeletonPartition(*g)); + BOOST_CHECK(out->is()); + BOOST_CHECK_EQUAL(out->as().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 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 out(algorithm::straightSkeletonPartition(*g)); + BOOST_CHECK(out->is()); + BOOST_CHECK_EQUAL(out->as().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 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 out(algorithm::straightSkeletonPartition(*g)); + BOOST_CHECK(out->is()); + 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 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 out(algorithm::straightSkeletonPartition(*g)); + BOOST_CHECK(out->is()); + 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 g(io::readWkt("POLYGON EMPTY")); + std::unique_ptr out(algorithm::straightSkeletonPartition(*g)); + BOOST_CHECK(out->is()); + 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 g(io::readWkt("LINESTRING(0 0, 1 1, 2 2)")); + BOOST_CHECK_THROW(algorithm::straightSkeletonPartition(*g), std::exception); +} + BOOST_AUTO_TEST_SUITE_END() diff --git a/test/unit/SFCGAL/capi/sfcgal_cTest.cpp b/test/unit/SFCGAL/capi/sfcgal_cTest.cpp index 43c37e2b..7a87d26d 100644 --- a/test/unit/SFCGAL/capi/sfcgal_cTest.cpp +++ b/test/unit/SFCGAL/capi/sfcgal_cTest.cpp @@ -237,4 +237,29 @@ BOOST_AUTO_TEST_CASE(testForceRHR_3D) BOOST_CHECK_EQUAL(expectedGeom, strApi); delete[] wkbApi; } + +BOOST_AUTO_TEST_CASE(testStraightSkeletonPartitionC) +{ + + sfcgal_set_error_handlers(printf, on_error); + std::unique_ptr 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()