diff --git a/CHANGELOG.md b/CHANGELOG.md index d0567cdf..7c9adf29 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,6 +28,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Removed * Remove unnecessary headers from `compas_cgal.h`. +* Remove `reconstruction.h` and `reconstruction.cpp`. +* Remove `straight_skeleton_2.h` and `straight_skeleton_2.cpp`. ## [0.7.3] 2025-03-18 diff --git a/CMakeLists.txt b/CMakeLists.txt index 07c5f725..3f0c302a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -13,6 +13,16 @@ set(CMAKE_CXX_FLAGS_DEBUG "") # ===================================================================== option(ENABLE_PRECOMPILED_HEADERS "Enable precompiled headers for the build" OFF) + +# ===================================================================== +# Set maximum heap size for MSVC +# ===================================================================== + +if(MSVC) +set(CMAKE_GENERATOR_PLATFORM x64) + add_compile_options(/Zm1200) +endif() + # ===================================================================== # Build size reduction. # ===================================================================== diff --git a/docs/api.rst b/docs/api.rst index 0a6b6c4f..33ac51df 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -9,10 +9,8 @@ API Reference api/compas_cgal.intersections api/compas_cgal.measure api/compas_cgal.meshing - 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 diff --git a/docs/api/compas_cgal.reconstruction.rst b/docs/api/compas_cgal.reconstruction.rst deleted file mode 100644 index 185f241f..00000000 --- a/docs/api/compas_cgal.reconstruction.rst +++ /dev/null @@ -1,15 +0,0 @@ -******************************************************************************** -compas_cgal.reconstruction -******************************************************************************** - -.. currentmodule:: compas_cgal.reconstruction - -.. autosummary:: - :toctree: generated/ - :nosignatures: - - pointset_normal_estimation - pointset_outlier_removal - pointset_reduction - pointset_smoothing - poisson_surface_reconstruction diff --git a/docs/api/compas_cgal.straight_skeleton_2.rst b/docs/api/compas_cgal.straight_skeleton_2.rst deleted file mode 100644 index 7daa1bfb..00000000 --- a/docs/api/compas_cgal.straight_skeleton_2.rst +++ /dev/null @@ -1,15 +0,0 @@ -******************************************************************************** -compas_cgal.straight_skeleton_2 -******************************************************************************** - -.. currentmodule:: compas_cgal.straight_skeleton_2 - -.. autosummary:: - :toctree: generated/ - :nosignatures: - - interior_straight_skeleton - interior_straight_skeleton_with_holes - offset_polygon - weighted_offset_polygon - offset_polygon_with_holes diff --git a/src/compas_cgal.cpp b/src/compas_cgal.cpp index 382fed94..f49609b4 100644 --- a/src/compas_cgal.cpp +++ b/src/compas_cgal.cpp @@ -6,12 +6,12 @@ void init_meshing(nb::module_ &); void init_measure(nb::module_ &); void init_booleans(nb::module_ &); void init_intersections(nb::module_ &); -void init_reconstruction(nb::module_ &); +// void init_reconstruction(nb::module_ &); void init_skeletonization(nb::module_ &); void init_slicer(nb::module_ &); void init_subdivision(nb::module_ &); void init_triangulation(nb::module_ &); -void init_straight_skeleton_2(nb::module_ &); +// void init_straight_skeleton_2(nb::module_ &); /** * @brief Module initialization function for COMPAS CGAL extension @@ -47,7 +47,7 @@ NB_MODULE(compas_cgal_ext, m) { nb::bind_vector>(m, "VectorRowMatrixXd"); // Be aware that both Pybind11 and Nanobind makes copy for vectors init_intersections(m); - init_reconstruction(m); +// init_reconstruction(m); init_skeletonization(m); @@ -59,5 +59,5 @@ NB_MODULE(compas_cgal_ext, m) { init_triangulation(m); - init_straight_skeleton_2(m); +// init_straight_skeleton_2(m); } diff --git a/src/compas_cgal/reconstruction.py b/src/compas_cgal/reconstruction.py deleted file mode 100644 index f79694dc..00000000 --- a/src/compas_cgal/reconstruction.py +++ /dev/null @@ -1,170 +0,0 @@ -from typing import Tuple -from typing import Union - -import numpy as np -from compas.geometry import Point -from compas.geometry import Vector - -from compas_cgal.compas_cgal_ext import reconstruction - -from .types import FloatNx3 -from .types import IntNx3 - - -def poisson_surface_reconstruction( - points: Union[list[Point], FloatNx3], - normals: Union[list[Vector], FloatNx3], -) -> Tuple[FloatNx3, IntNx3]: - """Reconstruct a surface from a point cloud using the Poisson surface reconstruction algorithm. - - Parameters - ---------- - points : list of :class:`compas.geometry.Point` or :class:`numpy.ndarray` - The points of the point cloud. - normals : list of :class:`compas.geometry.Vector` or :class:`numpy.ndarray` - The normals of the point cloud. - - Returns - ------- - tuple of :class:`numpy.ndarray` - The vertices and faces of the reconstructed surface. - - Raises - ------ - ValueError - If points or normals are not 3D - If number of points and normals don't match - If less than 3 points are provided - If points are not well-distributed for reconstruction - RuntimeError - If the reconstruction algorithm fails - - Notes - ----- - The Poisson surface reconstruction algorithm requires: - 1. A sufficiently dense point cloud - 2. Well-oriented normals - 3. Points distributed across a meaningful surface - """ - # Convert input to numpy arrays with proper type and memory layout - P = np.asarray(points, dtype=np.float64, order="C") - N = np.asarray(normals, dtype=np.float64, order="C") - - # Input validation - if P.ndim != 2 or P.shape[1] != 3: - raise ValueError("Points must be an Nx3 array") - if N.ndim != 2 or N.shape[1] != 3: - raise ValueError("Normals must be an Nx3 array") - if P.shape[0] != N.shape[0]: - raise ValueError("Number of points and normals must match") - if len(P) < 3: - raise ValueError("At least 3 points are required") - - # Ensure normals are unit vectors - norms = np.linalg.norm(N, axis=1) - if not np.allclose(norms, 1.0): - N = N / norms[:, np.newaxis] - - try: - return reconstruction.poisson_surface_reconstruction(P, N) - except RuntimeError as e: - raise RuntimeError(f"Poisson surface reconstruction failed: {str(e)}") - - -def pointset_outlier_removal( - points: Union[list[Point], FloatNx3], - nnnbrs: int = 10, - radius: float = 1.0, -) -> FloatNx3: - """Remove outliers from a point cloud using the point set outlier removal algorithm. - - Parameters - ---------- - points : list of :class:`compas.geometry.Point` or :class:`numpy.ndarray` - The points of the point cloud. - nnnbrs : int, optional - The number of nearest neighbors to consider for each point. - radius : float, optional - The radius of the sphere to consider for each point as a multiplication factor of the average point spacing. - - Returns - ------- - :class:`numpy.ndarray` - The points of the point cloud without outliers. - - """ - P = np.asarray(points, dtype=np.float64, order="C") - return reconstruction.pointset_outlier_removal(P, nnnbrs, radius) - - -def pointset_reduction( - points: Union[list[Point], FloatNx3], - spacing: float = 2, -) -> FloatNx3: - """Remove outliers from a point cloud using the point set outlier removal algorithm. - - Parameters - ---------- - points : list of :class:`compas.geometry.Point` or :class:`numpy.ndarray` - The points of the point cloud. - spacing : int, optional - The cell size. - - Returns - ------- - :class:`numpy.ndarray` - The vectors of the point cloud. - - """ - P = np.asarray(points, dtype=np.float64, order="C") - return reconstruction.pointset_reduction(P, spacing) - - -def pointset_smoothing( - points: Union[list[Point], FloatNx3], - neighbors: int = 8, - iterations: int = 1, -) -> FloatNx3: - """Remove outliers from a point cloud using the point set outlier removal algorithm. - - Parameters - ---------- - points : list of :class:`compas.geometry.Point` or :class:`numpy.ndarray` - The points of the point cloud. - neighbors : int, optional - The number of nearest neighbors to consider for each point. - - Returns - ------- - :class:`numpy.ndarray` - The vectors of the point cloud. - - """ - P = np.asarray(points, dtype=np.float64, order="C") - return reconstruction.pointset_smoothing(P, neighbors, iterations) - - -def pointset_normal_estimation( - points: Union[list[Point], FloatNx3], - neighbors: int = 8, - erase: bool = False, -) -> Tuple[FloatNx3, FloatNx3]: - """Remove outliers from a point cloud using the point set outlier removal algorithm. - - Parameters - ---------- - points : list of :class:`compas.geometry.Point` or :class:`numpy.ndarray` - The points of the point cloud. - neighbors : int, optional - The number of nearest neighbors to consider for each point. - erase : bool, optional - Erase points that are not oriented properly. - - Returns - ------- - :class:`numpy.ndarray` - The vectors of the point cloud. - - """ - P = np.asarray(points, dtype=np.float64, order="C") - return reconstruction.pointset_normal_estimation(P, neighbors, erase) diff --git a/src/compas_cgal/straight_skeleton_2.py b/src/compas_cgal/straight_skeleton_2.py deleted file mode 100644 index 04e11266..00000000 --- a/src/compas_cgal/straight_skeleton_2.py +++ /dev/null @@ -1,252 +0,0 @@ -from typing import Tuple -from typing import Union - -import numpy as np -from compas.datastructures import Graph -from compas.geometry import Polygon -from compas.geometry import normal_polygon -from compas.tolerance import TOL - -from compas_cgal.compas_cgal_ext import straight_skeleton_2 - -from .types import IntNx1 -from .types import IntNx2 -from .types import VerticesNumpy - - -def graph_from_skeleton_data(points: VerticesNumpy, indices: IntNx1, edges: IntNx2, edge_types: IntNx1) -> Graph: - """Create a graph from the skeleton data. - - Parameters - ---------- - points : :class:`numpy.ndarray` - The vertices of the skeleton, each vertex defined by 3 spatial coordinates. - indices : :class:`numpy.ndarray` - The vertex indices of the skeleton, corresponding to the points. - edges : :class:`numpy.ndarray` - The edges of the skeleton, each edge defined by 2 vertex indices. - edge_types : :class:`numpy.ndarray` - The type per edge, `0` for inner bisector, `1` for bisector, and `2` for boundary. - - Returns - ------- - :class:`compas.datastructures.Graph` - The skeleton as a graph. - """ - graph = Graph() - for pt, i in zip(points, indices): - graph.add_node(key=i, x=pt[0], y=pt[1], z=pt[2]) - - for edge, etype in zip(edges, edge_types): - edge = graph.add_edge(*edge) - if etype == 0: - graph.edge_attribute(edge, "inner_bisector", True) - elif etype == 1: - graph.edge_attribute(edge, "bisector", True) - else: - graph.edge_attribute(edge, "boundary", True) - return graph - - -def interior_straight_skeleton(points, as_graph=True) -> Union[Graph, Tuple[VerticesNumpy, IntNx1, IntNx2, IntNx1]]: - """Compute the skeleton of a 2D polygon. - - Parameters - ---------- - points : list of point coordinates or :class:`compas.geometry.Polygon` - The points of the polygon. - as_graph : bool, optional - Whether the skeleton should be returned as a graph, defaults to `True`. - - Returns - ------- - :attr:`compas.datastructures.Graph` or tuple of (vertices, indices, edges, edge_types) - The skeleton of the polygon. - - Raises - ------ - ValueError - If the normal of the polygon is not directed vertically upwards like [0, 0, 1]. - """ - points = list(points) - normal = normal_polygon(points, True) - if not TOL.is_allclose(normal, [0, 0, 1]): - raise ValueError("The normal of the polygon should be [0, 0, 1]. The normal of the provided polygon is {}".format(normal)) - V = np.asarray(points, dtype=np.float64, order="C") - points, indices, edges, edge_types = straight_skeleton_2.create_interior_straight_skeleton(V) - if as_graph: - return graph_from_skeleton_data(points, indices, edges, edge_types) - return points, indices, edges, edge_types - - -def interior_straight_skeleton_with_holes(points, holes, as_graph=True) -> Union[Graph, Tuple[VerticesNumpy, IntNx1, IntNx2, IntNx1]]: - """Compute the skeleton of a 2D polygon with holes. - - Parameters - ---------- - points : list of point coordinates or :class:`compas.geometry.Polygon` - The points of the 2D polygon. - holes : list of list of point coordinates or list of :class:`compas.geometry.Polygon` - The holes of the polygon. - as_graph : bool, optional - Whether the skeleton should be returned as a graph, defaults to `True`. - - Returns - ------- - :attr:`compas.datastructures.Graph` or tuple of (vertices, indices, edges, edge_types) - The skeleton of the polygon. - - Raises - ------ - ValueError - If the normal of the polygon is not directed vertically upwards like [0, 0, 1]. - If the normal of a hole is not directed vertically downwards like [0, 0, -1]. - """ - points = list(points) - normal = normal_polygon(points, True) - if not TOL.is_allclose(normal, [0, 0, 1]): - raise ValueError("The normal of the polygon should be [0, 0, 1]. The normal of the provided polygon is {}".format(normal)) - V = np.asarray(points, dtype=np.float64, order="C") - - H = [] - for i, hole in enumerate(holes): - points = list(hole) - normal_hole = normal_polygon(points, True) - if not TOL.is_allclose(normal_hole, [0, 0, -1]): - raise ValueError("The normal of the hole should be [0, 0, -1]. The normal of the provided {}-th hole is {}".format(i, normal_hole)) - hole = np.asarray(points, dtype=np.float64) - H.append(hole) - points, indices, edges, edge_types = straight_skeleton_2.create_interior_straight_skeleton_with_holes(V, H) - if as_graph: - return graph_from_skeleton_data(points, indices, edges, edge_types) - return points, indices, edges, edge_types - - -def offset_polygon(points, offset) -> list[Polygon]: - """Compute the offset from a 2D polygon. - - Parameters - ---------- - points : list of point coordinates or :class:`compas.geometry.Polygon` - The points of the 2D polygon. - offset : float - The offset distance. If negative, the offset is outside the polygon, otherwise inside. - - Returns - ------- - list[:class:`Polygon`] - The offset polygon(s). - - Raises - ------ - ValueError - If the normal of the polygon is not directed vertically upwards like [0, 0, 1]. - """ - points = list(points) - normal = normal_polygon(points, True) - if not TOL.is_allclose(normal, [0, 0, 1]): - raise ValueError("The normal of the polygon should be [0, 0, 1]. The normal of the provided polygon is {}".format(normal)) - V = np.asarray(points, dtype=np.float64, order="C") - offset = float(offset) - if offset < 0: # outside - offset_polygons = straight_skeleton_2.create_offset_polygons_2_outer(V, abs(offset))[1:] # first one is box - else: # inside - offset_polygons = straight_skeleton_2.create_offset_polygons_2_inner(V, offset) - return [Polygon(points.tolist()) for points in offset_polygons] - - -def offset_polygon_with_holes(points, holes, offset) -> list[Tuple[Polygon, list[Polygon]]]: - """Compute the offset from a 2D polygon with holes. - - Parameters - ---------- - points : list of point coordinates or :class:`compas.geometry.Polygon` - The points of the 2D polygon. - holes : list of list of point coordinates or list of :class:`compas.geometry.Polygon` - The holes of the polygon. - offset : float - The offset distance. If negative, the offset is outside the polygon, otherwise inside. - - Returns - ------- - list of tuple of (:class:`Polygon`, list[:class:`Polygon`]) - The polygons with holes. - - Raises - ------ - ValueError - If the normal of the polygon is not directed vertically upwards like [0, 0, 1]. - If the normal of a hole is not directed vertically downwards like [0, 0, -1]. - """ - points = list(points) - normal = normal_polygon(points, True) - if not TOL.is_allclose(normal, [0, 0, 1]): - raise ValueError("The normal of the polygon should be [0, 0, 1]. The normal of the provided polygon is {}".format(normal)) - V = np.asarray(points, dtype=np.float64, order="C") - - H = [] - for i, hole in enumerate(holes): - points = hole - normal_hole = normal_polygon(points, True) - if not TOL.is_allclose(normal_hole, [0, 0, -1]): - raise ValueError("The normal of the hole should be [0, 0, -1]. The normal of the provided {}-th hole is {}".format(i, normal_hole)) - hole = np.asarray(points, dtype=np.float64, order="C") - H.append(hole) - - if offset < 0: # outside - offset_polygons = straight_skeleton_2.create_offset_polygons_2_outer_with_holes(V, H, abs(offset)) - else: # inside - offset_polygons = straight_skeleton_2.create_offset_polygons_2_inner_with_holes(V, H, offset) - - result = [] - for points, holes_np in offset_polygons: - polygon = Polygon(points.tolist()) - holes = [] - for hole in holes_np: - holes.append(Polygon(hole.tolist())) - result.append((polygon, holes)) - return result - - -def weighted_offset_polygon(points, offset, weights) -> list[Polygon]: - """Compute the offset from a 2D polygon with weights. - - Parameters - ---------- - points : list of point coordinates or :class:`compas.geometry.Polygon` - The points of the 2D polygon. - offset : float - The offset distance. If negative, the offset is outside the polygon, otherwise inside. - weights : list of float - The weights for each edge, starting with the edge between the last and the first point. - - Returns - ------- - list[:class:`Polygon`] - The offset polygon(s). - - Raises - ------ - ValueError - If the normal of the polygon is not directed vertically upwards like [0, 0, 1]. - ValueError - If the number of weights does not match the number of points. - """ - points = list(points) - normal = normal_polygon(points, True) - if not TOL.is_allclose(normal, [0, 0, 1]): - raise ValueError("The normal of the polygon should be [0, 0, 1]. The normal of the provided polygon is {}".format(normal)) - - V = np.asarray(points, dtype=np.float64, order="C") - offset = float(offset) - # Reshape weights to be a column vector (n x 1) - W = np.asarray(weights, dtype=np.float64, order="C").reshape(-1, 1) - if W.shape[0] != V.shape[0]: - raise ValueError("The number of weights should be equal to the number of points %d != %d." % (W.shape[0], V.shape[0])) - if offset < 0: - offset_polygons = straight_skeleton_2.create_weighted_offset_polygons_2_outer(V, abs(offset), W) - # Return all except the first polygon which is the bounding box - return [Polygon(points.tolist()) for points in offset_polygons[1:]] - else: - offset_polygons = straight_skeleton_2.create_weighted_offset_polygons_2_inner(V, offset, W) - return [Polygon(points.tolist()) for points in offset_polygons] diff --git a/src/reconstruction.cpp b/src/reconstruction.cpp deleted file mode 100644 index ea223a93..00000000 --- a/src/reconstruction.cpp +++ /dev/null @@ -1,279 +0,0 @@ -#include "reconstruction.h" - -typedef std::pair PointVectorPair; -typedef CGAL::Point_set_3 PointSet; -typedef CGAL::Parallel_if_available_tag ConcurrencyTag; - -std::tuple -poisson_surface_reconstruction( - Eigen::Ref points, - Eigen::Ref normals) -{ - compas::Polyhedron mesh; - std::vector points_with_normals; - - for (int i = 0; i < points.rows(); i++) - { - points_with_normals.push_back(PointVectorPair( - compas::Point(points(i, 0), points(i, 1), points(i, 2)), - compas::Vector(normals(i, 0), normals(i, 1), normals(i, 2)))); - } - - double average_spacing = CGAL::compute_average_spacing( - points_with_normals, - 6, - CGAL::parameters::point_map(CGAL::First_of_pair_property_map())); - - CGAL::poisson_surface_reconstruction_delaunay( - points_with_normals.begin(), - points_with_normals.end(), - CGAL::First_of_pair_property_map(), - CGAL::Second_of_pair_property_map(), - mesh, - average_spacing); - - return compas::polyhedron_to_vertices_and_faces(mesh); -} - - -compas::RowMatrixXd -pointset_outlier_removal( - Eigen::Ref points, - int num_neighbors, - double radius) -{ - int p = points.rows(); - std::vector points_vector; - - for (int i = 0; i < p; i++) - { - points_vector.push_back(compas::Point(points(i, 0), points(i, 1), points(i, 2))); - } - - const double average_spacing = CGAL::compute_average_spacing(points_vector, num_neighbors); - - std::vector::iterator to_remove = CGAL::remove_outliers( - points_vector, - num_neighbors, - CGAL::parameters::threshold_percent(100.).threshold_distance(radius * average_spacing)); - - points_vector.erase(to_remove, points_vector.end()); - std::vector(points_vector).swap(points_vector); - - std::size_t s = points_vector.size(); - compas::RowMatrixXd result(s, 3); - - for (int i = 0; i < s; i++) - { - result(i, 0) = (double)points_vector[i].x(); - result(i, 1) = (double)points_vector[i].y(); - result(i, 2) = (double)points_vector[i].z(); - } - - return result; -} - - -compas::RowMatrixXd -pointset_reduction( - Eigen::Ref points, - double spacing - ) -{ - - // ====================================================================== - // Convert Eigen matrix to a vector of cgal points. - // ====================================================================== - int p = points.rows(); - std::vector points_vector; - - for (int i = 0; i < p; i++) - points_vector.emplace_back(points.row(i)[0], points.row(i)[1], points.row(i)[2]); - - // ====================================================================== - // Simplification by voxelization. - // ====================================================================== - double cell_size = spacing; - auto iterator_to_first_to_remove = CGAL::grid_simplify_point_set(points_vector, cell_size); - points_vector.erase(iterator_to_first_to_remove, points_vector.end()); - - // Optional: after erase(), shrink_to_fit to trim excess capacity - points_vector.shrink_to_fit(); - - // ====================================================================== - // Convert the vector of cgal points to an Eigen matrix. - // ====================================================================== - std::size_t s = points_vector.size(); - compas::RowMatrixXd smoothed_eigen_points(s, 3); - - for (std::size_t i = 0; i < s; i++) - smoothed_eigen_points.row(i) << static_cast(points_vector[i].x()), static_cast(points_vector[i].y()), static_cast(points_vector[i].z()); - - return smoothed_eigen_points; - - -} - - -compas::RowMatrixXd -pointset_smoothing( - Eigen::Ref points, - int num_neighbors, - int num_iterations) -{ - - // ====================================================================== - // Convert Eigen matrix to a vector of cgal points. - // ====================================================================== - int p = points.rows(); - std::vector points_vector; - - for (int i = 0; i < p; i++) - points_vector.emplace_back(points.row(i)[0], points.row(i)[1], points.row(i)[2]); - - // ====================================================================== - // Smooth the points, default is 24. - // ====================================================================== - for (int i = 0; i < num_iterations; i++) - CGAL::jet_smooth_point_set(points_vector, num_neighbors); - - // ====================================================================== - // Convert the vector of cgal points to an Eigen matrix - // ====================================================================== - std::size_t s = points_vector.size(); - compas::RowMatrixXd smoothed_eigen_points(s, 3); - - for (std::size_t i = 0; i < s; i++) - smoothed_eigen_points.row(i) << static_cast(points_vector[i].x()), static_cast(points_vector[i].y()), static_cast(points_vector[i].z()); - - return smoothed_eigen_points; - -} - -std::tuple -pointset_normal_estimation( - Eigen::Ref points, - int num_neighbors, - bool erase_unoriented) -{ - - // ====================================================================== - // Convert Eigen matrix to CGAL point vector pair. - // ====================================================================== - int p = points.rows(); - std::list points_with_normals; - - for (int i = 0; i < p; i++) - points_with_normals.emplace_back(PointVectorPair( - compas::Point(points(i, 0), points(i, 1), points(i, 2)), - compas::Vector(0, 0, 0) - )); - - - // =========================================================================== - // Estimates normals direction. - // Case 1 - the radius is precomputed by the average spacing. - // Case 2 - the normals are estimated with a fixed number of neighbors. - // =========================================================================== - if ( num_neighbors < 1 ) { - - // First compute a spacing using the K parameter - double spacing - = CGAL::compute_average_spacing - (points_with_normals, num_neighbors, - CGAL::parameters::point_map(CGAL::First_of_pair_property_map())); - - // Then, estimate normals with a fixed radius - CGAL::pca_estimate_normals (points_with_normals,0, // when using a neighborhood radius, K=0 means no limit on the number of neighbors returns - CGAL::parameters::point_map (CGAL::First_of_pair_property_map ()). - normal_map (CGAL::Second_of_pair_property_map ()). - neighbor_radius (2. * spacing)); // use 2*spacing as neighborhood radius - - } else { - - CGAL::pca_estimate_normals - (points_with_normals, num_neighbors, - CGAL::parameters::point_map (CGAL::First_of_pair_property_map ()). - normal_map (CGAL::Second_of_pair_property_map ())); - } - - // =========================================================================== - // Orients normals, give at least 3 points to make a plane. - // =========================================================================== - unsigned orientation_neighbors = (num_neighbors < 3) ? 3 : num_neighbors; - std::list::iterator unoriented_points_begin = - CGAL::mst_orient_normals (points_with_normals, num_neighbors, - CGAL::parameters::point_map (CGAL::First_of_pair_property_map ()). - normal_map (CGAL::Second_of_pair_property_map ())); - - // =========================================================================== - // Optional: delete points with an unoriented normal. - // =========================================================================== - if(erase_unoriented) - points_with_normals.erase (unoriented_points_begin, points_with_normals.end ()); - - // =========================================================================== - // Convert CGAL normals vector pair to Eigen matrix. - // =========================================================================== - - std::size_t s = points_with_normals.size(); - compas::RowMatrixXd eigen_points(s, 3); - compas::RowMatrixXd eigen_vectors(s, 3); - - std::size_t i = 0; - for (const auto& e : points_with_normals) { - eigen_points.row(i) << static_cast(e.first.x()), static_cast(e.first.y()), static_cast(e.first.z()); - eigen_vectors.row(i) << static_cast(e.second.x()), static_cast(e.second.y()), static_cast(e.second.z()); - i++; - } - - return std::make_tuple(eigen_points, eigen_vectors); - -} - -void init_reconstruction(nb::module_& m) { - auto submodule = m.def_submodule("reconstruction"); - - submodule.def( - "poisson_surface_reconstruction", - &poisson_surface_reconstruction, - "Perform Poisson surface reconstruction on an oriented pointcloud with normals", - "points"_a, - "normals"_a - ); - - submodule.def( - "pointset_outlier_removal", - &pointset_outlier_removal, - "Remove outliers from a pointcloud", - "points"_a, - "num_neighbors"_a, - "radius"_a - ); - - submodule.def( - "pointset_reduction", - &pointset_reduction, - "Reduce number of points using hierarchy simplification", - "points"_a, - "spacing"_a = 2.0 - ); - - submodule.def( - "pointset_smoothing", - &pointset_smoothing, - "Smooth a pointcloud using jet smoothing", - "points"_a, - "num_neighbors"_a = 8, - "num_iterations"_a = 1 - ); - - submodule.def( - "pointset_normal_estimation", - &pointset_normal_estimation, - "Estimate pointcloud normals and orient them", - "points"_a, - "num_neighbors"_a = 8, - "erase_unoriented"_a = true - ); -} diff --git a/src/reconstruction.h b/src/reconstruction.h deleted file mode 100644 index a7190601..00000000 --- a/src/reconstruction.h +++ /dev/null @@ -1,91 +0,0 @@ -#pragma once - -#include "compas.h" - -// CGAL reconstruction -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -/** - * @brief Perform Poisson surface reconstruction on an oriented pointcloud with normals. - * - * @param points Matrix of point positions as Nx3 matrix in row-major order (float64) - * @param normals Matrix of point normals as Nx3 matrix in row-major order (float64) - * @return std::tuple containing: - * - vertices as Rx3 matrix (float64) - * - faces as Sx3 matrix (int32) - */ -std::tuple -poisson_surface_reconstruction( - Eigen::Ref points, - Eigen::Ref normals); - -/** - * @brief Remove outliers from a pointcloud. - * - * @param points Matrix of point positions as Nx3 matrix in row-major order (float64) - * @param num_neighbors The number of nearest neighbors to consider - * @param radius The radius of the sphere to consider as a multiplication factor of the average spacing - * @return compas::RowMatrixXd - * - Filtered points as Mx3 matrix (float64), where M <= N - */ -compas::RowMatrixXd -pointset_outlier_removal( - Eigen::Ref points, - int num_neighbors, - double radius); - -/** - * @brief Reduce number of points using hierarchy simplification. - * https://doc.cgal.org/latest/Point_set_processing_3/Point_set_processing_3_2hierarchy_simplification_example_8cpp-example.html - * - * @param points Matrix of point positions as Nx3 matrix in row-major order (float64) - * @param spacing The scale for the average spacing - * @return compas::RowMatrixXd - * - Reduced points as Mx3 matrix (float64), where M <= N - */ -compas::RowMatrixXd -pointset_reduction( - Eigen::Ref points, - double spacing = 2.0); - -/** - * @brief Smooth a pointcloud by number of neighbors and iterations using jet_smoothing algorithm. - * https://doc.cgal.org/4.3/Point_set_processing_3/Point_set_processing_3_2jet_smoothing_example_8cpp-example.html - * - * @param points Matrix of point positions as Nx3 matrix in row-major order (float64) - * @param num_neighbors The number of nearest neighbors to consider - * @param num_iterations The number of iterations - * @return compas::RowMatrixXd - * - Smoothed points as Nx3 matrix (float64) - */ -compas::RowMatrixXd -pointset_smoothing( - Eigen::Ref points, - int num_neighbors = 8, - int num_iterations = 1); - -/** - * @brief Estimate pointcloud normals and oriented them. - * https://doc.cgal.org/latest/Point_set_processing_3/Point_set_processing_3_2normals_example_8cpp-example.html - * - * @param points Matrix of point positions as Nx3 matrix in row-major order (float64) - * @param num_neighbors The number of nearest neighbors to consider - * @param erase_unoriented Erase points that are not oriented properly - * @return std::tuple containing: - * - points as Mx3 matrix (float64), where M <= N if erase_unoriented is true - * - normals as Mx3 matrix (float64) - */ -std::tuple -pointset_normal_estimation( - Eigen::Ref points, - int num_neighbors=8, - bool erase_unoriented=true); diff --git a/src/straight_skeleton_2.cpp b/src/straight_skeleton_2.cpp deleted file mode 100644 index ccde7fe5..00000000 --- a/src/straight_skeleton_2.cpp +++ /dev/null @@ -1,449 +0,0 @@ -#include "straight_skeleton_2.h" - -typedef CGAL::Exact_predicates_inexact_constructions_kernel K; -typedef K::Point_2 Point; -typedef CGAL::Polygon_2 Polygon_2; -typedef CGAL::Polygon_with_holes_2 Polygon_with_holes; -typedef CGAL::Straight_skeleton_2 Ss; -typedef Ss::Halfedge_const_handle Halfedge_const_handle; -typedef Ss::Vertex_const_handle Vertex_const_handle; - -// Change from boost to std shared_ptr -typedef std::shared_ptr SsPtr; -typedef std::shared_ptr PolygonPtr; -typedef std::vector PolygonPtrVector; -typedef std::shared_ptr PolygonWithHolesPtr; -typedef std::vector PolygonWithHolesPtrVector; - -std::tuple, compas::RowMatrixXi, std::vector> -mesh_data_from_skeleton(std::shared_ptr& skeleton) -{ - std::size_t num_vertices = skeleton->size_of_vertices(); - std::size_t num_edges = skeleton->size_of_halfedges() / 2; // halfedges are stored twice - - compas::RowMatrixXd vertex_matrix(num_vertices, 3); - std::vector vertex_indices; // to save the vertex ids - compas::RowMatrixXi edge_matrix(num_edges, 2); - std::vector edge_types; // to save the edge type: 0: inner bisector, 1: bisector, 2: boundary - - std::size_t i = 0; - for (auto vertex_iter = skeleton->vertices_begin(); vertex_iter != skeleton->vertices_end(); ++vertex_iter) - { - const Vertex_const_handle vertex = vertex_iter; - vertex_matrix(i, 0) = (double)vertex->point().x(); - vertex_matrix(i, 1) = (double)vertex->point().y(); - vertex_matrix(i, 2) = 0; - vertex_indices.push_back((int)vertex->id()); - i++; - } - i = 0; - for (auto edge_iter = skeleton->halfedges_begin(); edge_iter != skeleton->halfedges_end(); ++edge_iter) - { - const Halfedge_const_handle halfedge = edge_iter; - const Vertex_const_handle& vertex1 = halfedge->vertex(); - const Vertex_const_handle& vertex2 = halfedge->opposite()->vertex(); - - if (&*vertex1 < &*vertex2) - { - edge_matrix(i, 0) = (int)vertex1->id(); - edge_matrix(i, 1) = (int)vertex2->id(); - - if (halfedge->is_inner_bisector()) - { - edge_types.push_back(0); - } - else if (halfedge->is_bisector()) - { - edge_types.push_back(1); - } - else - { - edge_types.push_back(2); - } - i++; - } - } - return std::make_tuple(vertex_matrix, vertex_indices, edge_matrix, edge_types); -} - -compas::RowMatrixXd -polygon_to_data(Polygon_2 const& polygon) -{ - std::size_t n = polygon.size(); - compas::RowMatrixXd points(n, 3); - int i = 0; - for(typename Polygon_2::Vertex_const_iterator vertex_iter = polygon.vertices_begin(); - vertex_iter != polygon.vertices_end(); ++vertex_iter) - { - points(i, 0) = (double)(*vertex_iter).x(); - points(i, 1) = (double)(*vertex_iter).y(); - points(i, 2) = 0; - i++; - } - return points; -} - -std::tuple> -polygon_with_holes_to_data(Polygon_with_holes const& polygon_with_holes) -{ - std::vector holes; - compas::RowMatrixXd boundary_points = polygon_to_data(polygon_with_holes.outer_boundary()); - for(typename Polygon_with_holes::Hole_const_iterator hole_iter = polygon_with_holes.holes_begin(); - hole_iter != polygon_with_holes.holes_end(); ++hole_iter) - { - compas::RowMatrixXd hole_points = polygon_to_data(*hole_iter); - holes.push_back(hole_points); - } - return std::make_tuple(boundary_points, holes); -} - -Polygon_2 -data_to_polygon(const compas::RowMatrixXd& vertices) -{ - Polygon_2 polygon; - for (int i = 0; i < vertices.rows(); i++) { - polygon.push_back(Point(vertices(i, 0), vertices(i, 1))); - } - return polygon; -} - -Polygon_with_holes -data_to_polygon_with_holes( - const compas::RowMatrixXd& boundary_vertices, - const std::vector& hole_vertices) -{ - Polygon_2 outer_polygon = data_to_polygon(boundary_vertices); - Polygon_with_holes polygon_with_holes(outer_polygon); - for (const auto& hole_points : hole_vertices) { - Polygon_2 hole; - for (int i = 0; i < hole_points.rows(); i++) { - hole.push_back(Point(hole_points(i, 0), hole_points(i, 1))); - } - polygon_with_holes.add_hole(hole); - } - return polygon_with_holes; -} - -std::tuple, compas::RowMatrixXi, std::vector> -pmp_create_interior_straight_skeleton( - const compas::RowMatrixXd& vertices) -{ - Polygon_2 polygon = data_to_polygon(vertices); - SsPtr skeleton = CGAL::create_interior_straight_skeleton_2(polygon.vertices_begin(), polygon.vertices_end()); - return mesh_data_from_skeleton(skeleton); -} - -std::tuple, compas::RowMatrixXi, std::vector> -pmp_create_interior_straight_skeleton_with_holes( - const compas::RowMatrixXd& boundary_vertices, - const std::vector& hole_vertices) -{ - Polygon_with_holes polygon = data_to_polygon_with_holes(boundary_vertices, hole_vertices); - SsPtr skeleton = CGAL::create_interior_straight_skeleton_2(polygon); - return mesh_data_from_skeleton(skeleton); -} - -std::vector -pmp_create_offset_polygons_2_inner( - const compas::RowMatrixXd& vertices, - double& offset_distance) -{ - Polygon_2 polygon = data_to_polygon(vertices); - PolygonPtrVector offset_polygons = CGAL::create_interior_skeleton_and_offset_polygons_2(offset_distance, polygon); - - std::vector result; - for(const auto& offset_polygon : offset_polygons) { - compas::RowMatrixXd points = polygon_to_data(*offset_polygon); - result.push_back(points); - } - return result; -} - -std::vector>> -pmp_create_offset_polygons_2_inner_with_holes( - const compas::RowMatrixXd& boundary_vertices, - const std::vector& hole_vertices, - double& offset_distance) -{ - Polygon_with_holes polygon = data_to_polygon_with_holes(boundary_vertices, hole_vertices); - PolygonWithHolesPtrVector offset_polygons = CGAL::create_interior_skeleton_and_offset_polygons_with_holes_2( - offset_distance, polygon); - - std::vector>> result; - for(const auto& offset_polygon : offset_polygons) { - result.push_back(polygon_with_holes_to_data(*offset_polygon)); - } - return result; -} - -std::vector -pmp_create_offset_polygons_2_outer( - const compas::RowMatrixXd& vertices, - double& offset_distance) -{ - Polygon_2 polygon = data_to_polygon(vertices); - PolygonPtrVector offset_polygons = CGAL::create_exterior_skeleton_and_offset_polygons_2(offset_distance, polygon); - - std::vector result; - for(const auto& offset_polygon : offset_polygons) { - compas::RowMatrixXd points = polygon_to_data(*offset_polygon); - result.push_back(points); - } - return result; -} - -std::vector>> -pmp_create_offset_polygons_2_outer_with_holes( - const compas::RowMatrixXd& boundary_vertices, - const std::vector& hole_vertices, - double& offset_distance) -{ - Polygon_with_holes polygon = data_to_polygon_with_holes(boundary_vertices, hole_vertices); - PolygonWithHolesPtrVector offset_polygons = CGAL::create_exterior_skeleton_and_offset_polygons_with_holes_2( - offset_distance, polygon); - - std::vector>> result; - for(const auto& offset_polygon : offset_polygons) { - result.push_back(polygon_with_holes_to_data(*offset_polygon)); - } - return result; -} - -std::vector -pmp_create_weighted_offset_polygons_2_inner( - const compas::RowMatrixXd& vertices, - double offset_distance, - const compas::RowMatrixXd& edge_weights) -{ - if (edge_weights.rows() != vertices.rows()) { - throw std::invalid_argument("Number of weights must match number of polygon vertices"); - } - - Polygon_2 polygon = data_to_polygon(vertices); - std::vector weight_values; - weight_values.reserve(edge_weights.rows()); - for (int i = 0; i < edge_weights.rows(); i++) { - if (edge_weights(i, 0) <= 0) { - throw std::invalid_argument("Weights must be positive"); - } - weight_values.push_back(edge_weights(i, 0)); - } - - try { - SsPtr skeleton = CGAL::create_interior_weighted_straight_skeleton_2(polygon, weight_values); - if (!skeleton) { - throw std::runtime_error("Failed to create weighted straight skeleton"); - } - - PolygonPtrVector offset_polygons = CGAL::create_offset_polygons_2(offset_distance, *skeleton); - std::vector result; - result.reserve(offset_polygons.size()); - for(const auto& offset_polygon : offset_polygons) { - result.push_back(polygon_to_data(*offset_polygon)); - } - return result; - } catch (const CGAL::Precondition_exception& e) { - throw std::runtime_error("CGAL precondition failed: Invalid input for weighted offset"); - } -} - -std::vector -pmp_create_weighted_offset_polygons_2_outer( - const compas::RowMatrixXd& vertices, - double offset_distance, - const compas::RowMatrixXd& edge_weights) -{ - if (edge_weights.rows() != vertices.rows()) { - throw std::invalid_argument("Number of weights must match number of polygon vertices"); - } - - Polygon_2 polygon = data_to_polygon(vertices); - std::vector weight_values; - weight_values.reserve(edge_weights.rows()); - for (int i = 0; i < edge_weights.rows(); i++) { - if (edge_weights(i, 0) <= 0) { - throw std::invalid_argument("Weights must be positive"); - } - weight_values.push_back(edge_weights(i, 0)); - } - - try { - SsPtr skeleton = CGAL::create_exterior_weighted_straight_skeleton_2(offset_distance, weight_values, polygon); - if (!skeleton) { - throw std::runtime_error("Failed to create weighted straight skeleton"); - } - - PolygonPtrVector offset_polygons = CGAL::create_offset_polygons_2(offset_distance, *skeleton); - std::vector result; - result.reserve(offset_polygons.size()); - for(const auto& offset_polygon : offset_polygons) { - result.push_back(polygon_to_data(*offset_polygon)); - } - return result; - } catch (const CGAL::Precondition_exception& e) { - throw std::runtime_error("CGAL precondition failed: Invalid input for weighted offset"); - } -} - -void init_straight_skeleton_2(nb::module_& m) { - auto submodule = m.def_submodule("straight_skeleton_2"); - - submodule.def( - "create_interior_straight_skeleton", - &pmp_create_interior_straight_skeleton, - "Create an interior straight skeleton from a polygon.\n\n" - "Parameters\n" - "----------\n" - "vertices : array-like\n" - " Matrix of vertex positions (Nx2, float64)\n" - "\n" - "Returns\n" - "-------\n" - "tuple\n" - " - Matrix of skeleton vertices (Mx2, float64)\n" - " - Vector of source vertex indices\n" - " - Matrix of skeleton edges (Kx2, int32)\n" - " - Vector of source edge indices", - "vertices"_a); - - submodule.def( - "create_interior_straight_skeleton_with_holes", - &pmp_create_interior_straight_skeleton_with_holes, - "Create an interior straight skeleton from a polygon with holes.\n\n" - "Parameters\n" - "----------\n" - "boundary_vertices : array-like\n" - " Matrix of boundary vertex positions (Nx2, float64)\n" - "hole_vertices : list\n" - " List of hole vertex matrices (each Mx2, float64)\n" - "\n" - "Returns\n" - "-------\n" - "tuple\n" - " - Matrix of skeleton vertices (Px2, float64)\n" - " - Vector of source vertex indices\n" - " - Matrix of skeleton edges (Qx2, int32)\n" - " - Vector of source edge indices", - "boundary_vertices"_a, - "hole_vertices"_a); - - submodule.def( - "create_offset_polygons_2_inner", - &pmp_create_offset_polygons_2_inner, - "Create inward offset polygons from a simple polygon.\n\n" - "Parameters\n" - "----------\n" - "vertices : array-like\n" - " Matrix of vertex positions (Nx2, float64)\n" - "offset_distance : float\n" - " Offset distance (positive for inward)\n" - "\n" - "Returns\n" - "-------\n" - "list\n" - " List of offset polygon matrices (each Mx2, float64)", - "vertices"_a, - "offset_distance"_a); - - submodule.def( - "create_offset_polygons_2_inner_with_holes", - &pmp_create_offset_polygons_2_inner_with_holes, - "Create inward offset polygons from a polygon with holes.\n\n" - "Parameters\n" - "----------\n" - "boundary_vertices : array-like\n" - " Matrix of boundary vertex positions (Nx2, float64)\n" - "hole_vertices : list\n" - " List of hole vertex matrices (each Mx2, float64)\n" - "offset_distance : float\n" - " Offset distance (positive for inward)\n" - "\n" - "Returns\n" - "-------\n" - "list\n" - " List of tuples (outer polygon, list of hole polygons)", - "boundary_vertices"_a, - "hole_vertices"_a, - "offset_distance"_a); - - submodule.def( - "create_offset_polygons_2_outer", - &pmp_create_offset_polygons_2_outer, - "Create outward offset polygons from a simple polygon.\n\n" - "Parameters\n" - "----------\n" - "vertices : array-like\n" - " Matrix of vertex positions (Nx2, float64)\n" - "offset_distance : float\n" - " Offset distance (positive for outward)\n" - "\n" - "Returns\n" - "-------\n" - "list\n" - " List of offset polygon matrices (each Mx2, float64)", - "vertices"_a, - "offset_distance"_a); - - submodule.def( - "create_offset_polygons_2_outer_with_holes", - &pmp_create_offset_polygons_2_outer_with_holes, - "Create outward offset polygons from a polygon with holes.\n\n" - "Parameters\n" - "----------\n" - "boundary_vertices : array-like\n" - " Matrix of boundary vertex positions (Nx2, float64)\n" - "hole_vertices : list\n" - " List of hole vertex matrices (each Mx2, float64)\n" - "offset_distance : float\n" - " Offset distance (positive for outward)\n" - "\n" - "Returns\n" - "-------\n" - "list\n" - " List of tuples (outer polygon, list of hole polygons)", - "boundary_vertices"_a, - "hole_vertices"_a, - "offset_distance"_a); - - submodule.def( - "create_weighted_offset_polygons_2_inner", - &pmp_create_weighted_offset_polygons_2_inner, - "Create inward weighted offset polygons from a simple polygon.\n\n" - "Parameters\n" - "----------\n" - "vertices : array-like\n" - " Matrix of vertex positions (Nx2, float64)\n" - "offset_distance : float\n" - " Offset distance (must be positive)\n" - "edge_weights : array-like\n" - " Matrix of edge weights (Nx1, float64, must be positive)\n" - "\n" - "Returns\n" - "-------\n" - "list\n" - " List of offset polygon matrices (each Mx2, float64)", - "vertices"_a, - "offset_distance"_a, - "edge_weights"_a); - - submodule.def( - "create_weighted_offset_polygons_2_outer", - &pmp_create_weighted_offset_polygons_2_outer, - "Create outward weighted offset polygons from a simple polygon.\n\n" - "Parameters\n" - "----------\n" - "vertices : array-like\n" - " Matrix of vertex positions (Nx2, float64)\n" - "offset_distance : float\n" - " Offset distance (must be positive)\n" - "edge_weights : array-like\n" - " Matrix of edge weights (Nx1, float64, must be positive)\n" - "\n" - "Returns\n" - "-------\n" - "list\n" - " List of offset polygon matrices (each Mx2, float64)", - "vertices"_a, - "offset_distance"_a, - "edge_weights"_a); -} \ No newline at end of file diff --git a/src/straight_skeleton_2.h b/src/straight_skeleton_2.h deleted file mode 100644 index 126c0d7f..00000000 --- a/src/straight_skeleton_2.h +++ /dev/null @@ -1,139 +0,0 @@ -#pragma once - -#include "compas.h" - -// CGAL straight skeleton 2 -#include -#include -#include -#include -#include -#include -#include - - -/** - * @brief Creates a straight skeleton from a simple polygon without holes. - * - * @param vertices Matrix of polygon vertices as Nx2 matrix in row-major order (float64) - * @return std::tuple, compas::RowMatrixXi, std::vector> containing: - * - Matrix of skeleton vertices (Mx2, float64) - * - Vector of source vertex indices from input polygon for each skeleton vertex - * - Matrix of skeleton edges as vertex pairs (Kx2, int32) - * - Vector of source edge indices from input polygon for each skeleton edge - */ -std::tuple, compas::RowMatrixXi, std::vector> -pmp_create_interior_straight_skeleton( - const compas::RowMatrixXd& vertices -); - -/** - * @brief Creates a straight skeleton from a polygon with holes. - * - * @param boundary_vertices Matrix of boundary polygon vertices as Nx2 matrix in row-major order (float64) - * @param hole_vertices Vector of hole polygons, each as Mx2 matrix in row-major order (float64) - * @return std::tuple, compas::RowMatrixXi, std::vector> containing: - * - Matrix of skeleton vertices (Px2, float64) - * - Vector of source vertex indices from input polygon for each skeleton vertex - * - Matrix of skeleton edges as vertex pairs (Qx2, int32) - * - Vector of source edge indices from input polygon for each skeleton edge - */ -std::tuple, compas::RowMatrixXi, std::vector> -pmp_create_interior_straight_skeleton_with_holes( - const compas::RowMatrixXd& boundary_vertices, - const std::vector& hole_vertices -); - -/** - * @brief Creates offset polygons from a simple polygon. - * - * @param vertices Matrix of polygon vertices as Nx2 matrix in row-major order (float64) - * @param offset_distance Offset distance (positive for inward, negative for outward) - * @return std::vector Vector of offset polygons, each as Mx2 matrix of vertices (float64) - */ -std::vector -pmp_create_offset_polygons_2_inner( - const compas::RowMatrixXd& vertices, - double& offset_distance -); - -/** - * @brief Creates offset polygons from a polygon with holes. - * - * @param boundary_vertices Matrix of boundary polygon vertices as Nx2 matrix in row-major order (float64) - * @param hole_vertices Vector of hole polygons, each as Mx2 matrix in row-major order (float64) - * @param offset_distance Offset distance (positive for inward, negative for outward) - * @return std::vector>> Vector containing: - * - For each offset: tuple of outer polygon (Px2, float64) and vector of inner polygons (each Qx2, float64) - */ -std::vector>> -pmp_create_offset_polygons_2_inner_with_holes( - const compas::RowMatrixXd& boundary_vertices, - const std::vector& hole_vertices, - double& offset_distance -); - -/** - * @brief Creates outward offset polygons from a simple polygon. - * - * @param vertices Matrix of polygon vertices as Nx2 matrix in row-major order (float64) - * @param offset_distance Offset distance (positive for inward, negative for outward) - * @return std::vector Vector of offset polygons, each as Mx2 matrix of vertices (float64) - */ -std::vector -pmp_create_offset_polygons_2_outer( - const compas::RowMatrixXd& vertices, - double& offset_distance -); - -/** - * @brief Creates outward offset polygons from a polygon with holes. - * - * @param boundary_vertices Matrix of boundary polygon vertices as Nx2 matrix in row-major order (float64) - * @param hole_vertices Vector of hole polygons, each as Mx2 matrix in row-major order (float64) - * @param offset_distance Offset distance (positive for inward, negative for outward) - * @return std::vector>> Vector containing: - * - For each offset: tuple of outer polygon (Px2, float64) and vector of inner polygons (each Qx2, float64) - */ -std::vector>> -pmp_create_offset_polygons_2_outer_with_holes( - const compas::RowMatrixXd& boundary_vertices, - const std::vector& hole_vertices, - double& offset_distance -); - -/** - * @brief Create weighted offset polygons for the interior of a polygon. - * - * @param vertices Matrix of vertex coordinates as Nx2 matrix in row-major order (float64) - * @param offset_distance Offset distance (must be positive) - * @param edge_weights Matrix of weights for each edge as Nx1 matrix (float64, must be positive) - * @return std::vector Vector of offset polygons, each as Mx2 matrix of vertices (float64) - * @throws std::invalid_argument if edge_weights size doesn't match vertices or if weights are not positive - * @throws std::runtime_error if CGAL fails to create the skeleton or offset - */ -std::vector -pmp_create_weighted_offset_polygons_2_inner( - const compas::RowMatrixXd& vertices, - double offset_distance, - const compas::RowMatrixXd& edge_weights -); - -/** - * @brief Create weighted offset polygons for the exterior of a polygon. - * - * @param vertices Matrix of vertex coordinates as Nx2 matrix in row-major order (float64) - * @param offset_distance Offset distance (must be positive) - * @param edge_weights Matrix of weights for each edge as Nx1 matrix (float64, must be positive) - * @return std::vector Vector of offset polygons, each as Mx2 matrix of vertices (float64) - * @throws std::invalid_argument if edge_weights size doesn't match vertices or if weights are not positive - * @throws std::runtime_error if CGAL fails to create the skeleton or offset - */ -std::vector -pmp_create_weighted_offset_polygons_2_outer( - const compas::RowMatrixXd& vertices, - double offset_distance, - const compas::RowMatrixXd& edge_weights -); - -void init_straight_skeleton_2(nb::module_& m); \ No newline at end of file diff --git a/tests/test_reconstruction_pointset_normal_estimation.py b/tests/test_reconstruction_pointset_normal_estimation.py deleted file mode 100644 index 01435956..00000000 --- a/tests/test_reconstruction_pointset_normal_estimation.py +++ /dev/null @@ -1,30 +0,0 @@ -from pathlib import Path - -from compas.geometry import Pointcloud -from compas.geometry import Line -from compas_cgal.reconstruction import pointset_normal_estimation -from compas_cgal.reconstruction import pointset_reduction - - -def reconstruction_pointset_normal_estimation(): - FILE = Path(__file__).parent.parent / "data" / "box.ply" - - cloud = Pointcloud.from_ply(FILE) - reduced_cloud = Pointcloud(pointset_reduction(cloud, 1000)) - points, vectors = pointset_normal_estimation(reduced_cloud, 16, True) - - lines = [] - line_scale = 1000 - - for p, v in zip(points, vectors): - line = Line( - [p[0], p[1], p[2]], - [ - p[0] + v[0] * line_scale, - p[1] + v[1] * line_scale, - p[2] + v[2] * line_scale, - ], - ) - lines.append(line) - - assert len(lines) == len(points) diff --git a/tests/test_reconstruction_pointset_outlier_removal.py b/tests/test_reconstruction_pointset_outlier_removal.py deleted file mode 100644 index 8540dd1e..00000000 --- a/tests/test_reconstruction_pointset_outlier_removal.py +++ /dev/null @@ -1,40 +0,0 @@ -from pathlib import Path - -from compas.geometry import Pointcloud -from compas_cgal.reconstruction import pointset_outlier_removal -import numpy as np -from line_profiler import profile - - -@profile -def reconstruction_pointset_outlier_removal(): - FILE = Path(__file__).parent.parent / "data" / "forked_branch_1.ply" - c1 = Pointcloud.from_ply(FILE) - points = pointset_outlier_removal(c1, 30, 2.0) - c2 = Pointcloud(points) - c3 = c1.difference(c2) - c2.name = "pointset_outlier_removal" - c3.name = "pointset_outlier_removal" - return c2, c3 - - -def test_pointset_outlier_removal(): - # Create a point cloud with an obvious outlier - points = np.array( - [ - [0, 0, 0], - [1, 0, 0], - [0, 1, 0], - [1, 1, 0], - [10, 10, 10], # Outlier - ], - dtype=np.float64, - ) - - # Test outlier removal - filtered_points = pointset_outlier_removal(points, nnnbrs=3, radius=1.0) - - # Basic checks - assert filtered_points.shape[0] == 4 # Should have removed one point - assert filtered_points.shape[1] == 3 # Each point should be 3D - assert not np.any(np.all(filtered_points == [10, 10, 10], axis=1)) # Outlier should be removed diff --git a/tests/test_reconstruction_pointset_reduction.py b/tests/test_reconstruction_pointset_reduction.py deleted file mode 100644 index 42642961..00000000 --- a/tests/test_reconstruction_pointset_reduction.py +++ /dev/null @@ -1,18 +0,0 @@ -import numpy as np -from compas_cgal.reconstruction import pointset_reduction - - -def test_pointset_reduction(): - # Create a dense point cloud - x = np.linspace(0, 1, 10) - y = np.linspace(0, 1, 10) - X, Y = np.meshgrid(x, y) - points = np.column_stack((X.ravel(), Y.ravel(), np.zeros_like(X.ravel()))) - - # Test point reduction - spacing = 0.3 # Should reduce points significantly - reduced_points = pointset_reduction(points, spacing=spacing) - - # Basic checks - assert reduced_points.shape[0] < points.shape[0] # Should have fewer points - assert reduced_points.shape[1] == 3 # Each point should be 3D diff --git a/tests/test_reconstruction_pointset_smoothing.py b/tests/test_reconstruction_pointset_smoothing.py deleted file mode 100644 index c145c3d9..00000000 --- a/tests/test_reconstruction_pointset_smoothing.py +++ /dev/null @@ -1,22 +0,0 @@ -import numpy as np -from compas_cgal.reconstruction import pointset_smoothing - - -def test_pointset_smoothing(): - # Create a noisy point cloud - points = np.array( - [ - [0, 0, 0], - [1.1, 0.1, 0.1], # Noisy point - [0, 1, 0], - [0.9, 0.9, -0.1], # Noisy point - ], - dtype=np.float64, - ) - - # Test smoothing - smoothed_points = pointset_smoothing(points, neighbors=3, iterations=2) - - # Basic checks - assert smoothed_points.shape == points.shape # Should preserve point count - assert smoothed_points.shape[1] == 3 # Each point should be 3D diff --git a/tests/test_reconstruction_poisson_surface_reconstruction.py b/tests/test_reconstruction_poisson_surface_reconstruction.py deleted file mode 100644 index fdd8b560..00000000 --- a/tests/test_reconstruction_poisson_surface_reconstruction.py +++ /dev/null @@ -1,35 +0,0 @@ -import math -from pathlib import Path - -from compas.datastructures import Mesh -from compas.geometry import Pointcloud -from compas.geometry import Rotation -from compas.geometry import Scale -from compas_cgal.reconstruction import poisson_surface_reconstruction - - -def test_reconstruction_poisson_surface_reconstruction(): - FILE = Path(__file__).parent.parent / "data" / "oni.xyz" - - points = [] - normals = [] - with open(FILE, "r") as f: - for line in f: - x, y, z, nx, ny, nz = line.strip().split() - points.append([float(x), float(y), float(z)]) - normals.append([float(nx), float(ny), float(nz)]) - - V, F = poisson_surface_reconstruction(points, normals) - mesh = Mesh.from_vertices_and_faces(V, F) - - R = Rotation.from_axis_and_angle([1, 0, 0], math.radians(90)) - S = Scale.from_factors([5, 5, 5]) - T = R * S - - c = Pointcloud(V) - c.transform(T) - mesh.transform(T) - - assert mesh.is_manifold() - assert mesh.number_of_vertices() > 0 - assert mesh.number_of_faces() > 0 diff --git a/tests/test_straight_skeleton_2_interior_straight_skeleton.py b/tests/test_straight_skeleton_2_interior_straight_skeleton.py deleted file mode 100644 index da5d2c12..00000000 --- a/tests/test_straight_skeleton_2_interior_straight_skeleton.py +++ /dev/null @@ -1,32 +0,0 @@ -from compas_cgal.straight_skeleton_2 import interior_straight_skeleton - - -def test_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), - ] - - graph = interior_straight_skeleton(points) - - # Basic validation - edges = list(graph.edges()) - vertices = list(graph.nodes()) - assert len(edges) > 0 - assert len(vertices) > 0 - - # Check that we have some inner bisectors - has_inner_bisector = False - for edge in edges: - if graph.edge_attribute(edge, "inner_bisector"): - has_inner_bisector = True - break - assert has_inner_bisector diff --git a/tests/test_straight_skeleton_2_interior_straight_skeleton_offset_polygon.py b/tests/test_straight_skeleton_2_interior_straight_skeleton_offset_polygon.py deleted file mode 100644 index c10769a2..00000000 --- a/tests/test_straight_skeleton_2_interior_straight_skeleton_offset_polygon.py +++ /dev/null @@ -1,33 +0,0 @@ -from compas.geometry import Polygon -from compas_cgal.straight_skeleton_2 import offset_polygon - - -def test_offset_polygon(): - 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), - ] - - offset = 1.5 - - # Test inner offset - offset_polygon_inner = offset_polygon(points, offset) - assert len(offset_polygon_inner) > 0 - for polygon in offset_polygon_inner: - assert isinstance(polygon, Polygon) - assert len(polygon.points) > 0 - - # Test outer offset - offset_polygon_outer = offset_polygon(points, -offset) - assert len(offset_polygon_outer) > 0 - for polygon in offset_polygon_outer: - assert isinstance(polygon, Polygon) - assert len(polygon.points) > 0 diff --git a/tests/test_straight_skeleton_2_interior_straight_skeleton_weighted_offset_polygons.py b/tests/test_straight_skeleton_2_interior_straight_skeleton_weighted_offset_polygons.py deleted file mode 100644 index 297319fe..00000000 --- a/tests/test_straight_skeleton_2_interior_straight_skeleton_weighted_offset_polygons.py +++ /dev/null @@ -1,23 +0,0 @@ -from compas_cgal.straight_skeleton_2 import weighted_offset_polygon - - -def test_weighted_offset_polygon(): - 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), - ] - - distances = [0.1, 0.3, 0.6, 0.1, 0.7, 0.5, 0.2, 0.4, 0.8, 0.2] - weights = [1.0 / d for d in distances] - offset = 1.0 - offset_polygons_outer = weighted_offset_polygon(points, -offset, weights) - - assert len(offset_polygons_outer) == 1 diff --git a/tests/test_straight_skeleton_2_interior_straight_skeleton_with_holes.py b/tests/test_straight_skeleton_2_interior_straight_skeleton_with_holes.py deleted file mode 100644 index bec7e9fa..00000000 --- a/tests/test_straight_skeleton_2_interior_straight_skeleton_with_holes.py +++ /dev/null @@ -1,47 +0,0 @@ -from compas.geometry import Polygon -from compas_cgal.straight_skeleton_2 import interior_straight_skeleton_with_holes - - -def test_interior_straight_skeleton_with_holes(): - 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), - ] - - holes = [ - [(0.42, 0.88, 0.0), (1.1, -1.0, 0.0), (-1.97, -0.93, 0.0), (-1.25, 1.82, 0.0)], - [(4.25, -0.64, 0.0), (2.9, -3.03, 0.0), (2.12, -2.16, 0.0), (2.89, -0.36, 0.0)], - [(10.6, 0.29, 0.0), (9.48, -1.54, 0.0), (5.48, -1.26, 0.0), (5.98, -0.04, 0.0)], - ] - - polygon = Polygon(points) - holes = [Polygon(hole) for hole in holes] - graph = interior_straight_skeleton_with_holes(polygon, holes) - - # Basic validation - edges = list(graph.edges()) - vertices = list(graph.nodes()) - assert len(edges) > 0 - assert len(vertices) > 0 - - # Check that we have both inner bisectors and boundary edges - has_inner_bisector = False - has_boundary = False - for edge in edges: - if graph.edge_attribute(edge, "inner_bisector"): - has_inner_bisector = True - elif not graph.edge_attribute(edge, "bisector"): - has_boundary = True - if has_inner_bisector and has_boundary: - break - - assert has_inner_bisector - assert has_boundary