ADD poinset method: reduction, smoothing, normal estimation

docs/examples/pointset_reduction.py

@@ -0,0 +1,37 @@
+from pathlib import Path
+from compas.geometry import Pointcloud, Translation
+from compas_view2.app import App
+from compas_cgal.reconstruction import pointset_reduction
+
+# Define the file path for the point cloud data
+FILE = Path(__file__).parent.parent.parent / "data" / "forked_branch_1.ply"
+
+# Load the original point cloud
+original_points = Pointcloud.from_ply(FILE)
+
+# Create a copy of the point cloud for processing
+cloud = Pointcloud.from_ply(FILE)
+
+# Translate the original point cloud
+cloud.transform(Translation.from_vector([-1000, 0, 0]))
+
+# Apply point set reduction to the translated point cloud
+points = pointset_reduction(cloud, 50)
+print(f"Original points: {len(cloud)}, Reduced points: {len(points)}")
+
+# Initialize the COMPAS viewer
+viewer = App(width=1600, height=900)
+
+# Adjust viewer settings
+viewer.view.camera.scale = 1000
+viewer.view.grid.cell_size = 1000
+
+# Add the reduced point cloud and the original point cloud to the viewer
+viewer.add(Pointcloud(points))
+viewer.add(Pointcloud(original_points))
+
+# Set the camera to zoom to fit all points
+viewer.view.camera.zoom_extents()
+
+# Run the viewer
+viewer.run()

docs/examples/pointset_reduction.rst

@@ -0,0 +1,11 @@
+********************************************************************************
+Pointset Reduction
+********************************************************************************
+
+.. figure:: /_images/cgal_pointset_reduction.png
+    :figclass: figure
+    :class: figure-img img-fluid
+
+
+.. literalinclude:: pointset_reduction.py
+    :language: python

docs/examples/pointsets_normals.py

@@ -0,0 +1,51 @@
+from pathlib import Path
+from compas.geometry import Pointcloud, Point, Line
+from compas_view2.app import App
+from compas_view2.collections import Collection
+from compas.colors import Color
+from compas_cgal.reconstruction import pointset_normal_estimation
+
+# Define the file path for the point cloud data
+FILE = Path(__file__).parent.parent.parent / "data" / "forked_branch_1.ply"
+
+# Load the point cloud data from the PLY file
+cloud = Pointcloud.from_ply(FILE)
+
+# Estimate normals for the point cloud
+points, vectors = pointset_normal_estimation(cloud, 16, True)
+print(f"Original points: {len(cloud)}, Points with normals: {len(points)}, Vectors: {len(vectors)}")
+
+# Create lines and properties for visualizing normals
+lines = []
+line_properties = []
+line_scale = 25
+
+# Iterate over points and vectors to create lines and color properties
+for p, v in zip(points, vectors):
+    lines.append(Line(Point(p[0], p[1], p[2]), Point(p[0] + v[0] * line_scale, p[1] + v[1] * line_scale, p[2] + v[2] * line_scale)))
+
+    # Normalize vector components to be in the range [0, 1] for color representation
+    r = (v[0] + 1) * 0.5
+    g = (v[1] + 1) * 0.5
+    b = (v[2] + 1) * 0.5
+
+    # Store line color properties
+    line_properties.append({'linecolor': Color(r, g, b)})
+
+# Initialize the COMPAS viewer
+viewer = App(width=1600, height=900)
+
+# Adjust viewer settings
+viewer.view.camera.scale = 1000
+viewer.view.grid.cell_size = 1000
+
+# Create a line collection and add it to the viewer
+line_collection = Collection(lines, line_properties)
+viewer.add(Pointcloud(points))
+viewer.add(line_collection)
+
+# Set the camera to zoom to fit all points and lines
+viewer.view.camera.zoom_extents()
+
+# Run the viewer
+viewer.run()

docs/examples/pointsets_normals.rst

@@ -0,0 +1,11 @@
+********************************************************************************
+Pointset Normal Estimation
+********************************************************************************
+
+.. figure:: /_images/cgal_pointset_normal_estimation.png
+    :figclass: figure
+    :class: figure-img img-fluid
+
+
+.. literalinclude:: pointsets_normals.py
+    :language: python

docs/examples/pointsets_smoothing.py

@@ -0,0 +1,37 @@
+from pathlib import Path
+from compas.geometry import Pointcloud, Translation
+from compas_view2.app import App
+from compas_cgal.reconstruction import pointset_smoothing
+
+# Define the path to the PLY file
+ply_file_path = Path(__file__).parent.parent.parent / "data" / "box.ply"
+
+# Load the original point cloud and translate it
+original_points = Pointcloud.from_ply(ply_file_path)
+original_points.transform(Translation.from_vector([-10000, 0, 0]))
+
+# Load another copy of the point cloud for comparison and translate it in the opposite direction
+transformed_points = Pointcloud.from_ply(ply_file_path)
+transformed_points.transform(Translation.from_vector([10000, 0, 0]))
+
+# Apply point set smoothing to the transformed point cloud
+smoothed_points = pointset_smoothing(transformed_points, 1000, 3)
+
+# Create Pointcloud objects for visualization
+cloud_original = Pointcloud(original_points)
+cloud_transformed = Pointcloud(smoothed_points)
+
+# Set up the viewer
+viewer = App(width=1600, height=900)
+viewer.view.camera.scale = 1000
+viewer.view.grid.cell_size = 1000
+
+# Add the Pointclouds to the viewer
+viewer.add(cloud_original)
+viewer.add(cloud_transformed)
+
+# Adjust the camera and grid settings
+viewer.view.camera.zoom_extents()
+
+# Run the viewer
+viewer.run()

docs/examples/pointsets_smoothing.rst

@@ -0,0 +1,11 @@
+********************************************************************************
+Pointset Smoothing
+********************************************************************************
+
+.. figure:: /_images/cgal_pointset_smoothing.png
+    :figclass: figure
+    :class: figure-img img-fluid
+
+
+.. literalinclude:: pointsets_smoothing.py
+    :language: python

src/compas_cgal/reconstruction.py

@@ -55,3 +55,76 @@ def pointset_outlier_removal(
     """
     P = np.asarray(points, dtype=np.float64)
     return reconstruction.pointset_outlier_removal(P, nnnbrs, radius)
+
+
+def pointset_reduction(
+    points: Union[list[Point], NDArray[Shape["Any, 3"], Float]],
+    spacing: float = 2,
+) -> NDArray[Shape["Any, 3"], Float]:
+    """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)
+    return reconstruction.pointset_reduction(P, spacing)
+
+
+def pointset_smoothing(
+    points: Union[list[Point], NDArray[Shape["Any, 3"], Float]],
+    neighbors: int = 8,
+    iterations: int = 1,
+) -> NDArray[Shape["Any, 3"], Float]:
+    """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)
+    return reconstruction.pointset_smoothing(P, neighbors, iterations)
+
+
+def pointset_normal_estimation(
+    points: Union[list[Point], NDArray[Shape["Any, 3"], Float]],
+    neighbors: int = 8,
+    erase: bool = False,
+) -> Tuple[NDArray[Shape["Any, 3"], Float], NDArray[Shape["Any, 3"], Float]]:
+    """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)
+    return reconstruction.pointset_normal_estimation(P, neighbors, erase)