Added straight skeleton 2

CHANGELOG.md

@@ -9,6 +9,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+* Added `compas_cgal.straight_skeleton_2.create_interior_straight_skeleton`.
+
 ### Changed
 
 ### Removed

docs/api.rst

@@ -12,6 +12,7 @@ API Reference
     api/compas_cgal.reconstruction
     api/compas_cgal.slicer
     api/compas_cgal.skeletonization
+    api/compas_cgal.straight_skeleton_2
     api/compas_cgal.subdivision
     api/compas_cgal.triangulation
     api/compas_cgal.types

docs/api/compas_cgal.straight_skeleton_2.rst

@@ -0,0 +1,11 @@
+********************************************************************************
+compas_cgal.straight_skeleton_2
+********************************************************************************
+
+.. currentmodule:: compas_cgal.straight_skeleton_2
+
+.. autosummary::
+    :toctree: generated/
+    :nosignatures:
+
+    create_interior_straight_skeleton

docs/examples/straight_skeleton_2.py

@@ -0,0 +1,31 @@
+from compas.datastructures import Graph
+from compas.geometry import Polygon
+from compas_viewer import Viewer
+
+from compas_cgal.straight_skeleton_2 import create_interior_straight_skeleton
+
+points = [
+    (-1.91, 3.59, 0.0),
+    (-5.53, -5.22, 0.0),
+    (-0.39, -1.98, 0.0),
+    (2.98, -5.51, 0.0),
+    (4.83, -2.02, 0.0),
+    (9.70, -3.63, 0.0),
+    (12.23, 1.25, 0.0),
+    (3.42, 0.66, 0.0),
+    (2.92, 4.03, 0.0),
+    (-1.91, 3.59, 0.0),
+]
+polygon = Polygon(points)
+lines = create_interior_straight_skeleton(points)
+graph = Graph.from_lines(lines)
+
+# ==============================================================================
+# Viz
+# ==============================================================================
+
+viewer = Viewer(width=1600, height=900)
+viewer.renderer_config.show_grid = False
+viewer.scene.add(graph, edgecolor=(1.0, 0.0, 0.0))
+viewer.scene.add(polygon)
+viewer.show()

docs/examples/straight_skeleton_2.rst

@@ -0,0 +1,11 @@
+********************************************************************************
+2D Straight Skeleton
+********************************************************************************
+
+.. figure:: /_images/cgal_straight_skeleton_2.png
+    :figclass: figure
+    :class: figure-img img-fluid
+
+
+.. literalinclude:: straight_skeleton_2.py
+    :language: python

include/compas.h

@@ -24,6 +24,8 @@ namespace compas
     using Mesh = CGAL::Surface_mesh<Point>;
     using RowMatrixXd = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
     using RowMatrixXi = Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+    using Edge = std::tuple<std::vector<double>, std::vector<double>>;
+    using Edges = std::list<Edge>;
 
     Polyhedron polyhedron_from_vertices_and_faces(const RowMatrixXd &V, const RowMatrixXi &F);
 

setup.py

@@ -52,6 +52,7 @@ def get_scip_library():
                 "src/skeletonization.cpp",
                 "src/reconstruction.cpp",
                 "src/polygonal_surface_reconstruction.cpp",
+                "src/straight_skeleton_2.cpp",
             ]
         ),
         include_dirs=["./include", get_eigen_include(), get_pybind_include()],

src/compas_cgal.cpp

@@ -12,6 +12,7 @@ void init_subdivision(pybind11::module &);
 void init_skeletonization(pybind11::module &);
 void init_reconstruction(pybind11::module &);
 void init_polygonal_surface_reconstruction(pybind11::module &);
+void init_straight_skeleton_2(pybind11::module &);
 
 // the first parameter here ("_cgal") will be the name of the "so" or "pyd" file that will be produced by PyBind11
 // which is the entry point from where all other modules will be accessible.
@@ -30,4 +31,5 @@ PYBIND11_MODULE(_cgal, m)
     init_skeletonization(m);
     init_reconstruction(m);
     init_polygonal_surface_reconstruction(m);
+    init_straight_skeleton_2(m);
 }

src/compas_cgal/__init__.py

@@ -22,6 +22,7 @@
     "compas_cgal.measure",
     "compas_cgal.slicer",
     "compas_cgal.triangulation",
+    "compas_cgal.straight_skeleton_2",
 ]
 
 __all__ = ["HOME", "DATA", "DOCS", "TEMP"]

src/compas_cgal/straight_skeleton_2.py

@@ -0,0 +1,32 @@
+import numpy as np
+from compas.geometry import normal_polygon
+from compas.tolerance import TOL
+
+from compas_cgal._cgal import straight_skeleton_2
+
+from .types import PolylinesNumpy
+
+
+def create_interior_straight_skeleton(points) -> PolylinesNumpy:
+    """Compute the skeleton of a polygon.
+
+    Parameters
+    ----------
+    points : list of point coordinates or :class:`compas.geometry.Polygon`
+        The points of the polygon.
+
+    Returns
+    -------
+    :attr:`compas_cgal.types.PolylinesNumpy`
+        The skeleton of the polygon.
+
+    Raises
+    ------
+    ValueError
+        If the normal of the polygon is not [0, 0, 1].
+    """
+    points = list(points)
+    if not TOL.is_allclose(normal_polygon(points, True), [0, 0, 1]):
+        raise ValueError("Please pass a polygon with a normal vector of [0, 0, 1].")
+    V = np.asarray(points, dtype=np.float64)
+    return straight_skeleton_2.create_interior_straight_skeleton(V)

src/skeletonization.cpp

@@ -39,7 +39,7 @@ struct Split_polylines
     }
 };
 
-Edges pmp_mesh_skeleton(
+compas::Edges pmp_mesh_skeleton(
     Eigen::Ref<const compas::RowMatrixXd> &V,
     Eigen::Ref<const compas::RowMatrixXi> &F)
 {
@@ -63,7 +63,7 @@ Edges pmp_mesh_skeleton(
     // Split_polylines splitter(skeleton, polylines);
     // CGAL::split_graph_into_polylines(skeleton, splitter);
 
-    Edges edgelist;
+    compas::Edges edgelist;
 
     for (Skeleton_edge e : CGAL::make_range(edges(skeleton)))
     {
@@ -72,7 +72,7 @@ Edges pmp_mesh_skeleton(
 
         std::vector<double> s_vec = {s.x(), s.y(), s.z()};
         std::vector<double> t_vec = {t.x(), t.y(), t.z()};
-        Edge edge = std::make_tuple(s_vec, t_vec);
+        compas::Edge edge = std::make_tuple(s_vec, t_vec);
 
         edgelist.push_back(edge);
     }

src/skeletonization.h

@@ -3,10 +3,7 @@
 
 #include <compas.h>
 
-typedef std::tuple<std::vector<double>, std::vector<double>> Edge;
-typedef std::list<Edge> Edges;
-
-Edges pmp_mesh_skeleton(
+compas::Edges pmp_mesh_skeleton(
     Eigen::Ref<const compas::RowMatrixXd> &V,
     Eigen::Ref<const compas::RowMatrixXi> &F);
 

src/straight_skeleton_2.cpp

@@ -0,0 +1,55 @@
+
+#include "straight_skeleton_2.h"
+#include <CGAL/Polygon_2.h>
+#include <CGAL/create_straight_skeleton_2.h>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef K::Point_2 Point;
+typedef CGAL::Polygon_2<K> Polygon_2;
+typedef CGAL::Straight_skeleton_2<K> Ss;
+typedef boost::shared_ptr<Ss> SsPtr;
+typedef CGAL::Straight_skeleton_2<K>::Halfedge_const_handle Halfedge_const_handle;
+typedef CGAL::Straight_skeleton_2<K>::Vertex_const_handle Vertex_const_handle;
+
+compas::Edges pmp_create_interior_straight_skeleton(
+    Eigen::Ref<const compas::RowMatrixXd> &V)
+{
+    Polygon_2 poly;
+    for (int i = 0; i < V.rows(); i++)
+    {
+        poly.push_back(Point(V(i, 0), V(i, 1)));
+    }
+    SsPtr iss = CGAL::create_interior_straight_skeleton_2(poly.vertices_begin(), poly.vertices_end());
+    compas::Edges edgelist;
+    for(auto hit = iss->halfedges_begin(); hit != iss->halfedges_end(); ++hit){
+        const Halfedge_const_handle h = hit;
+        if(!h->is_bisector()){
+            continue;
+        }
+        const Vertex_const_handle& v1 = h->vertex();
+        const Vertex_const_handle& v2 = h->opposite()->vertex();
+        if(&*v1 < &*v2){
+            std::vector<double> s_vec = {v1->point().x(), v1->point().y(), 0};
+            std::vector<double> t_vec = {v2->point().x(), v2->point().y(), 0};
+            compas::Edge edge = std::make_tuple(s_vec, t_vec);
+            edgelist.push_back(edge);
+        }
+
+    }
+    return edgelist;
+};
+
+
+// ===========================================================================
+// PyBind11
+// ===========================================================================
+
+void init_straight_skeleton_2(pybind11::module &m)
+{
+    pybind11::module submodule = m.def_submodule("straight_skeleton_2");
+
+    submodule.def(
+        "create_interior_straight_skeleton",
+        &pmp_create_interior_straight_skeleton,
+        pybind11::arg("V").noconvert());
+};

src/straight_skeleton_2.h

@@ -0,0 +1,10 @@
+#ifndef COMPAS_STRAIGHT_SKELETON_2_H
+#define COMPAS_STRAIGHT_SKELETON_2_H
+
+#include <compas.h>
+
+
+compas::Edges pmp_create_interior_straight_skeleton(
+    Eigen::Ref<const compas::RowMatrixXd> &V);
+
+#endif /* COMPAS_STRAIGHT_SKELETON_2_H */

tests/test_straight_skeleton_2.py

@@ -0,0 +1,57 @@
+from compas_cgal.straight_skeleton_2 import create_interior_straight_skeleton
+from compas.tolerance import TOL
+
+
+def test_straight_polygon():
+    points = [
+        (-1, -1, 0),
+        (0, -12, 0),
+        (1, -1, 0),
+        (12, 0, 0),
+        (1, 1, 0),
+        (0, 12, 0),
+        (-1, 1, 0),
+        (-12, 0, 0),
+    ]
+    lines = create_interior_straight_skeleton(points)
+    assert len(lines) == 8
+
+
+def test_straight_polygon_2_compare():
+    points = [
+        (-1.91, 3.59, 0.0),
+        (-5.53, -5.22, 0.0),
+        (-0.39, -1.98, 0.0),
+        (2.98, -5.51, 0.0),
+        (4.83, -2.02, 0.0),
+        (9.70, -3.63, 0.0),
+        (12.23, 1.25, 0.0),
+        (3.42, 0.66, 0.0),
+        (2.92, 4.03, 0.0),
+        (-1.91, 3.59, 0.0),
+    ]
+    lines = create_interior_straight_skeleton(points)
+
+    expected = [
+        [[-1.91, 3.59, 0.0], [-0.139446292, 1.191439787, 0.0]],
+        [[-5.53, -5.22, 0.0], [-0.139446292, 1.191439787, 0.0]],
+        [[-0.39, -1.98, 0.0], [0.008499564, 1.241560466, 0.0]],
+        [[2.98, -5.51, 0.0], [2.44972507, -1.674799065, 0.0]],
+        [[4.83, -2.02, 0.0], [4.228131167, -0.522007766, 0.0]],
+        [[8.663865218, -1.084821998, 0.0], [9.7, -3.63, 0.0]],
+        [[12.23, 1.25, 0.0], [8.663865218, -1.084821998, 0.0]],
+        [[3.42, 0.66, 0.0], [1.755862468, -1.404991433, 0.0]],
+        [[2.92, 4.03, 0.0], [0.563706846, 1.033296141, 0.0]],
+        [[4.228131167, -0.522007766, 0.0], [2.44972507, -1.674799065, 0.0]],
+        [[4.228131167, -0.522007766, 0.0], [8.663865218, -1.084821998, 0.0]],
+        [[1.755862468, -1.404991433, 0.0], [2.44972507, -1.674799065, 0.0]],
+        [[0.563706846, 1.033296141, 0.0], [1.755862468, -1.404991433, 0.0]],
+        [[-0.139446292, 1.191439787, 0.0], [0.008499564, 1.241560466, 0.0]],
+        [[0.563706846, 1.033296141, 0.0], [0.008499564, 1.241560466, 0.0]],
+    ]
+    for act, exp in zip(lines, expected):
+        sa, ea = act
+        se, ee = exp
+        # the line direction sometimes changes ...
+        assert TOL.is_allclose(sa, se) or TOL.is_allclose(sa, ee)
+        assert TOL.is_allclose(ea, ee) or TOL.is_allclose(ea, se)