Extrude straight skeleton

Oslandia · Oct 10, 2023 · 47c7cbd · 47c7cbd
1 parent 35f72d2
commit 47c7cbd
Show file tree

Hide file tree

Showing 11 changed files with 537 additions and 146 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -68,7 +68,7 @@ if( ${CGAL_USE_AUTOLINK} )
 	add_definitions( "-DCGAL_NO_AUTOLINK" )
 endif()
 
-find_package( CGAL 5.3 COMPONENTS Core REQUIRED )
+find_package( CGAL 5.6 COMPONENTS Core REQUIRED )
 message( STATUS "CGAL ${CGAL_VERSION} found" )
 
 include_directories( ${CMAKE_BINARY_DIR}/include )

diff --git a/NEWS b/NEWS
@@ -1,6 +1,7 @@
 1.5.0 (2023-xx-xx):
     * Add polygon partition (Raphaël Delhome, Loïc Bartoletti)
     * WKT: Fix triangle code (Loïc Bartoletti)
+    * Straight Skeleton: Add a version with extrusion
 1.4.1 (2022-01-27):
     * Add alpha-shapes algorithm (Loïc Bartoletti, Hugo Mercier)
     * Fix build and tests for MSYS2/MinGW (Loïc Bartoletti)

diff --git a/src/PolyhedralSurface.cpp b/src/PolyhedralSurface.cpp
@@ -54,6 +54,24 @@ PolyhedralSurface::PolyhedralSurface(const MarkedPolyhedron &poly) : Surface()
   }
 }
 
+///
+///
+///
+PolyhedralSurface::PolyhedralSurface(const Mesh &sm) : Surface()
+{
+
+  using vertex_descriptor = Mesh::Vertex_index;
+  for (auto face : sm.faces()) {
+    auto *new_face = new LineString();
+    for (vertex_descriptor vd : vertices_around_face(sm.halfedge(face), sm)) {
+      new_face->addPoint(Point(sm.point(vd)));
+    }
+
+    new_face->addPoint(new_face->startPoint().clone());
+    _polygons.push_back(new Polygon(new_face));
+  }
+}
+
 ///
 ///
 ///

diff --git a/src/PolyhedralSurface.h b/src/PolyhedralSurface.h
@@ -11,12 +11,16 @@
 #include <boost/serialization/base_object.hpp>
 #include <vector>
 
+#include <SFCGAL/Kernel.h>
 #include <SFCGAL/Point.h>
 #include <SFCGAL/Polygon.h>
 #include <SFCGAL/TriangulatedSurface.h>
 #include <SFCGAL/triangulate/triangulatePolygon.h>
 
 #include <CGAL/Polyhedron_3.h>
+#include <CGAL/Surface_mesh.h>
+
+using Mesh = CGAL::Surface_mesh<SFCGAL::Kernel::Point_3>;
 
 namespace SFCGAL {
 
@@ -42,6 +46,10 @@ class SFCGAL_API PolyhedralSurface : public Surface {
    * Constructor from a CGAL::Polyhedron_3
    */
   PolyhedralSurface(const detail::MarkedPolyhedron &poly);
+  /**
+   * Constructor from a CGAL::Surface_mesh
+   */
+  PolyhedralSurface(const Mesh &sm);
   /**
    * Copy constructor
    */

diff --git a/src/algorithm/extrude.cpp b/src/algorithm/extrude.cpp
@@ -38,7 +38,8 @@ extrude(const Point &g, const Kernel::Vector_3 &v) -> LineString *;
 auto
 extrude(const LineString &g, const Kernel::Vector_3 &v) -> PolyhedralSurface *;
 auto
-extrude(const Polygon &g, const Kernel::Vector_3 &v) -> Solid *;
+extrude(const Polygon &g, const Kernel::Vector_3 &v, bool addTop = true)
+    -> Solid *;
 auto
 extrude(const Triangle &g, const Kernel::Vector_3 &v) -> Solid *;
 
@@ -120,7 +121,7 @@ extrude(const LineString &g, const Kernel::Vector_3 &v) -> PolyhedralSurface *
 ///
 ///
 auto
-extrude(const Polygon &g, const Kernel::Vector_3 &v) -> Solid *
+extrude(const Polygon &g, const Kernel::Vector_3 &v, bool addTop) -> Solid *
 {
   if (g.isEmpty()) {
     return new Solid();
@@ -142,11 +143,12 @@ extrude(const Polygon &g, const Kernel::Vector_3 &v) -> Solid *
   polyhedralSurface.addPolygon(bottom);
 
   // "top"
-  Polygon top(bottom);
-  top.reverse();
-  translate(top, v);
-  polyhedralSurface.addPolygon(top);
-
+  if (addTop) {
+    Polygon top(bottom);
+    top.reverse();
+    translate(top, v);
+    polyhedralSurface.addPolygon(top);
+  }
   // exterior ring and interior rings extruded
   for (size_t i = 0; i < bottom.numRings(); i++) {
     std::unique_ptr<PolyhedralSurface> boundaryExtruded(
@@ -391,5 +393,17 @@ extrude(const Geometry &g, const double &dx, const double &dy, const double &dz)
   return extrude(g, Kernel::FT(dx), Kernel::FT(dy), Kernel::FT(dz));
 }
 
+SFCGAL_API auto
+extrude(const Polygon &g, const double &height) -> std::unique_ptr<Geometry>
+{
+
+  if (!std::isfinite(height)) {
+    BOOST_THROW_EXCEPTION(NonFiniteValueException(
+        "trying to extrude with non finite value in direction"));
+  }
+
+  return std::unique_ptr<Geometry>(
+      extrude(g, Kernel::Vector_3(0.0, 0.0, height), false));
+}
 } // namespace algorithm
 } // namespace SFCGAL
diff --git a/src/algorithm/extrude.h b/src/algorithm/extrude.h
@@ -111,6 +111,8 @@ SFCGAL_API std::unique_ptr<Geometry>
 SFCGAL_API std::unique_ptr<Geometry>
            extrude(const Geometry &g, const Kernel::Vector_3 &v);
 
+SFCGAL_API std::unique_ptr<Geometry>
+           extrude(const Polygon &g, const double &height);
 } // namespace algorithm
 } // namespace SFCGAL
 

diff --git a/src/algorithm/straightSkeleton.cpp b/src/algorithm/straightSkeleton.cpp
@@ -10,27 +10,34 @@
 #include <SFCGAL/MultiLineString.h>
 #include <SFCGAL/MultiPolygon.h>
 #include <SFCGAL/Polygon.h>
+#include <SFCGAL/PolyhedralSurface.h>
+#include <SFCGAL/Solid.h>
 #include <SFCGAL/Triangle.h>
 
 #include <SFCGAL/Exception.h>
 
 #include <SFCGAL/algorithm/intersection.h>
 #include <SFCGAL/algorithm/isValid.h>
 #include <SFCGAL/algorithm/orientation.h>
+#include <SFCGAL/algorithm/tesselate.h>
 #include <SFCGAL/algorithm/translate.h>
 
 #include <CGAL/Straight_skeleton_converter_2.h>
+#include <CGAL/Surface_mesh.h>
 #include <CGAL/create_straight_skeleton_from_polygon_with_holes_2.h>
+#include <CGAL/extrude_skeleton.h>
 
 #include <memory>
 
 namespace SFCGAL {
 namespace algorithm {
 
 using Point_2              = Kernel::Point_2;
+using Point_3              = Kernel::Point_3;
 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>;
 
 namespace { // anonymous
 
@@ -379,5 +386,47 @@ approximateMedialAxis(const Geometry &g) -> std::unique_ptr<MultiLineString>
   return mx;
 }
 
+auto
+extrudeStraightSkeleton(const Polygon &g, double height)
+    -> std::unique_ptr<PolyhedralSurface>
+{
+
+  Mesh sm;
+  CGAL::extrude_skeleton(g.toPolygon_with_holes_2(), sm,
+                         CGAL::parameters::maximum_height(height));
+  std::unique_ptr<PolyhedralSurface> polys(new PolyhedralSurface(sm));
+
+  return polys;
+}
+
+auto
+extrudeStraightSkeleton(const Geometry &g, double height)
+    -> std::unique_ptr<PolyhedralSurface>
+{
+  SFCGAL_ASSERT_GEOMETRY_VALIDITY_2D(g);
+
+  if (g.geometryTypeId() != TYPE_POLYGON) {
+    BOOST_THROW_EXCEPTION(Exception("Geometry must be a Polygon"));
+  }
+  std::unique_ptr<PolyhedralSurface> result(
+      extrudeStraightSkeleton(g.as<Polygon>(), height));
+  propagateValidityFlag(*result, true);
+  return result;
+}
+
+auto
+extrudeStraightSkeleton(const Geometry &g, double building_height,
+                        double roof_height)
+    -> std::unique_ptr<PolyhedralSurface>
+{
+  std::unique_ptr<PolyhedralSurface> roof{
+      extrudeStraightSkeleton(g, roof_height)};
+  translate(*roof, 0.0, 0.0, building_height);
+  std::unique_ptr<Geometry> building(extrude(g.as<Polygon>(), building_height));
+  std::unique_ptr<PolyhedralSurface> result{
+      new PolyhedralSurface(building->as<Solid>().exteriorShell())};
+  result->addPolygons(*roof);
+  return result;
+};
 } // namespace algorithm
 } // namespace SFCGAL
diff --git a/src/algorithm/straightSkeleton.h b/src/algorithm/straightSkeleton.h
@@ -5,8 +5,10 @@
 #ifndef _SFCGAL_ALGORITHM_STRAIGHTSKELETON_H_
 #define _SFCGAL_ALGORITHM_STRAIGHTSKELETON_H_
 
+#include <SFCGAL/PolyhedralSurface.h>
+#include <SFCGAL/algorithm/extrude.h>
+#include <SFCGAL/algorithm/union.h>
 #include <SFCGAL/config.h>
-
 #include <memory>
 
 namespace SFCGAL {
@@ -84,6 +86,24 @@ SFCGAL_API std::unique_ptr<MultiLineString>
                             bool innerOnly = false, bool outputDistanceInM = false,
                             const double &toleranceAbs = 1e-8);
 
+/**
+ * @brief build a 3D straight skeleton extruded for a Polygon
+ * @ingroup detail
+ * @throws NotImplementedException If g is a Polygon with point touching rings.
+ */
+SFCGAL_API auto
+extrudedStraightSkeleton(const Polygon &g, double height)
+    -> std::unique_ptr<PolyhedralSurface>;
+
+SFCGAL_API auto
+extrudeStraightSkeleton(const Geometry &g, double height)
+    -> std::unique_ptr<PolyhedralSurface>;
+
+SFCGAL_API auto
+extrudeStraightSkeleton(const Geometry &g, double building_height,
+                        double roof_height)
+    -> std::unique_ptr<PolyhedralSurface>;
+
 } // namespace algorithm
 } // namespace SFCGAL
 

diff --git a/src/capi/sfcgal_c.cpp b/src/capi/sfcgal_c.cpp
@@ -1170,6 +1170,49 @@ sfcgal_geometry_has_validity_flag(const sfcgal_geometry_t *geom) -> int
   return g1->hasValidityFlag() ? 1 : 0;
 }
 
+extern "C" auto
+sfcgal_geometry_extrude_straight_skeleton(const sfcgal_geometry_t *geom,
+                                          double height) -> sfcgal_geometry_t *
+{
+  const auto *g1 = reinterpret_cast<const SFCGAL::Geometry *>(geom);
+  std::unique_ptr<SFCGAL::PolyhedralSurface> polys;
+
+  try {
+    polys = SFCGAL::algorithm::extrudeStraightSkeleton(*g1, height);
+  } catch (std::exception &e) {
+    SFCGAL_WARNING("During straight_extrude_skeleton_distance(A):");
+    SFCGAL_WARNING("  with A: %s",
+                   ((const SFCGAL::Geometry *)(geom))->asText().c_str());
+    SFCGAL_ERROR("%s", e.what());
+    return nullptr;
+  }
+
+  return polys.release();
+}
+
+extern "C" auto
+sfcgal_geometry_extrude_polygon_straight_skeleton(const sfcgal_geometry_t *geom,
+                                                  double building_height,
+                                                  double roof_height)
+    -> sfcgal_geometry_t *
+{
+  const auto *g1 = reinterpret_cast<const SFCGAL::Geometry *>(geom);
+  std::unique_ptr<SFCGAL::Geometry> polys;
+
+  try {
+    polys = SFCGAL::algorithm::extrudeStraightSkeleton(*g1, building_height,
+                                                       roof_height);
+  } catch (std::exception &e) {
+    SFCGAL_WARNING("During straight_extrude_skeleton_distance(A):");
+    SFCGAL_WARNING("  with A: %s",
+                   ((const SFCGAL::Geometry *)(geom))->asText().c_str());
+    SFCGAL_ERROR("%s", e.what());
+    return nullptr;
+  }
+
+  return polys.release();
+}
+
 extern "C" auto
 sfcgal_geometry_straight_skeleton_distance_in_m(const sfcgal_geometry_t *geom)
     -> sfcgal_geometry_t *
@@ -1259,15 +1302,15 @@ sfcgal_geometry_optimal_alpha_shapes(const sfcgal_geometry_t *geom,
   return result.release();
 }
 
-
 extern "C" sfcgal_geometry_t *
 sfcgal_y_monotone_partition_2(const sfcgal_geometry_t *geom)
 {
   const SFCGAL::Geometry *g1 = reinterpret_cast<const SFCGAL::Geometry *>(geom);
   std::unique_ptr<SFCGAL::Geometry> result;
 
   try {
-    result = SFCGAL::algorithm::partition_2(g1->as<const SFCGAL::Geometry>(), SFCGAL::algorithm::y_monotone);
+    result = SFCGAL::algorithm::partition_2(g1->as<const SFCGAL::Geometry>(),
+                                            SFCGAL::algorithm::y_monotone);
   } catch (std::exception &e) {
     SFCGAL_WARNING("During y_monotone_partition_2(A):");
     SFCGAL_WARNING("  with A: %s",
@@ -1286,7 +1329,8 @@ sfcgal_approx_convex_partition_2(const sfcgal_geometry_t *geom)
   std::unique_ptr<SFCGAL::Geometry> result;
 
   try {
-    result = SFCGAL::algorithm::partition_2(g1->as<const SFCGAL::Geometry>(), SFCGAL::algorithm::approx_convex);
+    result = SFCGAL::algorithm::partition_2(g1->as<const SFCGAL::Geometry>(),
+                                            SFCGAL::algorithm::approx_convex);
   } catch (std::exception &e) {
     SFCGAL_WARNING("During approx_convex_partition_2(A):");
     SFCGAL_WARNING("  with A: %s",
@@ -1305,7 +1349,9 @@ sfcgal_greene_approx_convex_partition_2(const sfcgal_geometry_t *geom)
   std::unique_ptr<SFCGAL::Geometry> result;
 
   try {
-    result = SFCGAL::algorithm::partition_2(g1->as<const SFCGAL::Geometry>(), SFCGAL::algorithm::greene_approx_convex);
+    result =
+        SFCGAL::algorithm::partition_2(g1->as<const SFCGAL::Geometry>(),
+                                       SFCGAL::algorithm::greene_approx_convex);
   } catch (std::exception &e) {
     SFCGAL_WARNING("During greene_approx_convex_partition_2(A):");
     SFCGAL_WARNING("  with A: %s",
@@ -1323,7 +1369,8 @@ sfcgal_optimal_convex_partition_2(const sfcgal_geometry_t *geom)
   std::unique_ptr<SFCGAL::Geometry> result;
 
   try {
-    result = SFCGAL::algorithm::partition_2(g1->as<const SFCGAL::Geometry>(), SFCGAL::algorithm::optimal_convex);
+    result = SFCGAL::algorithm::partition_2(g1->as<const SFCGAL::Geometry>(),
+                                            SFCGAL::algorithm::optimal_convex);
   } catch (std::exception &e) {
     SFCGAL_WARNING("During optimal_convex_partition_2(A):");
     SFCGAL_WARNING("  with A: %s",

diff --git a/src/capi/sfcgal_c.h b/src/capi/sfcgal_c.h
@@ -57,7 +57,7 @@ typedef enum {
   //     TYPE_SURFACE             = 14, //abstract
   SFCGAL_TYPE_POLYHEDRALSURFACE   = 15,
   SFCGAL_TYPE_TRIANGULATEDSURFACE = 16,
-  SFCGAL_TYPE_TRIANGLE   = 17,
+  SFCGAL_TYPE_TRIANGLE            = 17,
 
   //-- not official codes
   SFCGAL_TYPE_SOLID      = 101,
@@ -980,6 +980,21 @@ sfcgal_geometry_straight_skeleton(const sfcgal_geometry_t *geom);
 SFCGAL_API sfcgal_geometry_t *
 sfcgal_geometry_straight_skeleton_distance_in_m(const sfcgal_geometry_t *geom);
 
+/**
+ * Returns the extrude straight skeleton of the given Polygon
+ * @pre geom must be a Polygon
+ * @pre isValid(geom) == true
+ * @post isValid(return) == true
+ * @ingroup capi
+ */
+SFCGAL_API sfcgal_geometry_t *
+sfcgal_geometry_extrude_straight_skeleton(const sfcgal_geometry_t *geom,
+                                          double                   height);
+SFCGAL_API sfcgal_geometry_t *
+sfcgal_geometry_extrude_polygon_straight_skeleton(const sfcgal_geometry_t *geom,
+                                                  double building_height,
+                                                  double roof_height);
+
 /**
  * Returns the approximate medial axis for the given Polygon
  * Approximate medial axis is based on straight skeleton
@@ -1064,7 +1079,8 @@ SFCGAL_API sfcgal_geometry_t *
 sfcgal_approx_convex_partition_2(const sfcgal_geometry_t *geom);
 
 /**
- * Returns the greene approximal convex partition of a geometry (polygon without hole)
+ * Returns the greene approximal convex partition of a geometry (polygon without
+ * hole)
  * @pre isValid(geom) == true
  * @post isValid(return) == true
  * @ingroup capi