ADD intersections

Author: petrasvestartas

CMakeLists.txt

@@ -227,3 +227,4 @@ endfunction()
 # Add new modules here
 add_nanobind_module(booleans_ext src/booleans.cpp)
 add_nanobind_module(meshing_ext src/meshing.cpp)
+add_nanobind_module(intersections_ext src/intersections.cpp)

docs/examples/example_intersections.py

@@ -0,0 +1,74 @@
+import profile
+
+from compas.datastructures import Mesh
+from compas.geometry import Box
+from compas.geometry import Point
+from compas.geometry import Polyline
+from compas.geometry import Sphere
+from compas_viewer import Viewer
+from line_profiler import profile
+
+from compas_cgal.intersections import intersection_mesh_mesh
+
+
+@profile
+def main():
+    # ==============================================================================
+    # Make a box and a sphere
+    # ==============================================================================
+
+    box = Box(2)
+    A = box.to_vertices_and_faces(triangulated=True)
+
+    sphere = Sphere(1, point=[1, 1, 1])
+    B = sphere.to_vertices_and_faces(u=32, v=32, triangulated=True)
+
+    # ==============================================================================
+    # Compute the intersections
+    # ==============================================================================
+
+    pointsets = intersection_mesh_mesh(A, B)
+
+    # ==============================================================================
+    # Process output
+    # ==============================================================================
+
+    polylines = []
+    for points in pointsets:
+        points = [Point(*point) for point in points]
+        polyline = Polyline(points)
+        polylines.append(polyline)
+
+    return A, B, polylines
+
+
+A, B, polylines = main()
+
+# =============================================================================
+# Visualize
+# ==============================================================================
+
+viewer = Viewer()
+
+viewer.renderer.camera.target = [0, 0, 0]
+viewer.renderer.camera.position = [4, -6, 3]
+
+viewer.scene.add(Mesh.from_vertices_and_faces(*A), facecolor=(1.0, 0.0, 0.0), show_points=False)
+viewer.scene.add(
+    Mesh.from_vertices_and_faces(*B),
+    facecolor=(0.0, 1.0, 0.0),
+    show_points=False,
+    opacity=0.3,
+)
+
+for polyline in polylines:
+    viewer.scene.add(
+        polyline,
+        linecolor=(0.0, 0.0, 1.0),
+        lineswidth=3,
+        pointcolor=(0.0, 0.0, 1.0),
+        pointsize=20,
+        show_points=True,
+    )
+
+viewer.show()

docs/examples/example_intersections.rst

@@ -0,0 +1,17 @@
+Mesh Intersections
+=====================
+
+This example demonstrates how to compute intersections between two triangle meshes using COMPAS CGAL.
+
+Key Features:
+
+* Computing intersection curves between meshes
+* Visualizing intersection results
+* Handling multiple intersection curves
+
+.. figure:: /_images/example_intersections.png
+    :figclass: figure
+    :class: figure-img img-fluid
+
+.. literalinclude:: example_intersections.py
+   :language: python

src/booleans.cpp

@@ -104,43 +104,3 @@ NB_MODULE(booleans_ext, m) {
         "VB"_a,
         "FB"_a);
 }
-
-void init_booleans(nb::module_& m) {
-    auto submodule = m.def_submodule("booleans");
-
-    submodule.def(
-        "boolean_union",
-        &pmp_boolean_union,
-        "Boolean Union from triangular mesh vertices and faces.",
-        "VA"_a,
-        "FA"_a,
-        "VB"_a,
-        "FB"_a);
-
-    submodule.def(
-        "boolean_difference",
-        &pmp_boolean_difference,
-        "Boolean Difference from triangular mesh vertices and faces.",
-        "VA"_a,
-        "FA"_a,
-        "VB"_a,
-        "FB"_a);
-
-    submodule.def(
-        "boolean_intersection",
-        &pmp_boolean_intersection,
-        "Boolean Intersection from triangular mesh vertices and faces.",
-        "VA"_a,
-        "FA"_a,
-        "VB"_a,
-        "FB"_a);
-
-    submodule.def(
-        "split",
-        &pmp_split,
-        "Boolean Split from triangular mesh vertices and faces.",
-        "VA"_a,
-        "FA"_a,
-        "VB"_a,
-        "FB"_a);
-}

src/compas_cgal/intersections.py

@@ -0,0 +1,53 @@
+import numpy as np
+from compas.plugins import plugin
+
+from compas_cgal import intersections_ext
+
+from .types import PolylinesNumpy
+from .types import VerticesFaces
+
+
+@plugin(category="intersections")
+def intersection_mesh_mesh(
+    A: VerticesFaces,
+    B: VerticesFaces,
+) -> PolylinesNumpy:
+    """Compute the intersection of tow meshes.
+
+    Parameters
+    ----------
+    A : :attr:`compas_cgal.types.VerticesFaces`
+        Mesh A.
+    B : :attr:`compas_cgal.types.VerticesFaces`
+        Mesh B.
+
+    Returns
+    -------
+    :attr:`compas_cgal.types.PolylinesNumpy`
+        A list of intersection polylines, with each polyline an array of points.
+
+    Examples
+    --------
+    >>> from compas.geometry import Box, Sphere, Polyline
+    >>> from compas_cgal.intersections import intersection_mesh_mesh
+
+    >>> box = Box(1)
+    >>> sphere = Sphere(0.5, point=[1, 1, 1])
+
+    >>> A = box.to_vertices_and_faces(triangulated=True)
+    >>> B = sphere.to_vertices_and_faces(u=32, v=32, triangulated=True)
+
+    >>> result = intersection_mesh_mesh(A, B)
+    >>> polylines = [Polyline(points) for points in result]
+
+    """
+    VA, FA = A
+    VB, FB = B
+    VA = np.asarray(VA, dtype=np.float64, order="C")
+    FA = np.asarray(FA, dtype=np.int32, order="C")
+    VB = np.asarray(VB, dtype=np.float64, order="C")
+    FB = np.asarray(FB, dtype=np.int32, order="C")
+
+    pointsets: intersections_ext.VectorRowMatrixXd = intersections_ext.intersection_mesh_mesh(VA, FA, VB, FB)
+
+    return pointsets

src/intersections.cpp

@@ -0,0 +1,32 @@
+#include "intersections.h"
+
+std::vector<compas::RowMatrixXd>
+pmp_intersection_mesh_mesh(
+    Eigen::Ref<const compas::RowMatrixXd> vertices_a,
+    Eigen::Ref<const compas::RowMatrixXi> faces_a,
+    Eigen::Ref<const compas::RowMatrixXd> vertices_b,
+    Eigen::Ref<const compas::RowMatrixXi> faces_b)
+{
+    compas::Mesh mesh_a = compas::mesh_from_vertices_and_faces(vertices_a, faces_a);
+    compas::Mesh mesh_b = compas::mesh_from_vertices_and_faces(vertices_b, faces_b);
+
+    compas::Polylines polylines;
+    CGAL::Polygon_mesh_processing::surface_intersection(mesh_a, mesh_b, std::back_inserter(polylines));
+
+    std::vector<compas::RowMatrixXd> result = compas::polylines_to_lists_of_points(polylines);
+    return result;
+}
+
+NB_MODULE(intersections_ext, m) {
+
+    nb::bind_vector<std::vector<compas::RowMatrixXd>>(m, "VectorRowMatrixXd"); // Be aware that both Pybind11 and Nanobind makes copy for vectors
+
+    m.def(
+        "intersection_mesh_mesh",
+        &pmp_intersection_mesh_mesh,
+        "Compute intersection polylines between two triangle meshes.",
+        "vertices_a"_a,
+        "faces_a"_a,
+        "vertices_b"_a,
+        "faces_b"_a);
+}

src/intersections.h

@@ -0,0 +1,24 @@
+#pragma once
+
+#include "compas.h"
+
+// CGAL intersection
+#include <CGAL/Polygon_mesh_processing/intersection.h>
+
+/**
+ * Compute intersection between two triangle meshes.
+ *
+ * @param vertices_a Vertices of mesh A as Nx3 matrix in row-major order (float64)
+ * @param faces_a Faces of mesh A as Mx3 matrix of vertex indices in row-major order (int32)
+ * @param vertices_b Vertices of mesh B as Px3 matrix in row-major order (float64)
+ * @param faces_b Faces of mesh B as Qx3 matrix of vertex indices in row-major order (int32)
+ * @return std::vector<RowMatrixXd> containing:
+ *         - List of intersection polylines, each as Rx3 matrix of points (float64)
+ * @note Input meshes must be manifold and closed
+ */
+std::vector<compas::RowMatrixXd>
+pmp_intersection_mesh_mesh(
+    Eigen::Ref<const compas::RowMatrixXd> vertices_a,
+    Eigen::Ref<const compas::RowMatrixXi> faces_a,
+    Eigen::Ref<const compas::RowMatrixXd> vertices_b,
+    Eigen::Ref<const compas::RowMatrixXi> faces_b);

tests/test_intersections.py

@@ -0,0 +1,37 @@
+import pytest
+import numpy as np
+
+from compas.geometry import Box, Sphere
+from compas_cgal.intersections import intersection_mesh_mesh
+
+from compas_cgal.intersections_ext import VectorRowMatrixXd
+
+
+@pytest.fixture
+def box_sphere_meshes():
+    """Create test meshes: a box and a sphere that intersect."""
+    box = Box(2)  # 2x2x2 box centered at origin
+    box_mesh = box.to_vertices_and_faces(triangulated=True)
+
+    sphere = Sphere(1, point=[1, 1, 1])  # Sphere radius 1, centered at [1,1,1]
+    sphere_mesh = sphere.to_vertices_and_faces(u=32, v=32, triangulated=True)
+
+    return box_mesh, sphere_mesh
+
+
+def test_intersection_mesh_mesh(box_sphere_meshes):
+    """Test mesh-mesh intersection computation."""
+    mesh_a, mesh_b = box_sphere_meshes
+
+    # Compute intersections
+    pointsets = intersection_mesh_mesh(mesh_a, mesh_b)
+
+    # Basic validation
+    assert isinstance(pointsets, VectorRowMatrixXd), "Result should be a nanobind type"
+    assert len(pointsets) > 0, "Should find at least one intersection curve"
+
+    # Validate each pointset
+    for points in pointsets:
+        assert isinstance(points, np.ndarray), "Each pointset should be a numpy array"
+        assert points.shape[1] == 3, "Points should be 3D"
+        assert len(points) > 1, "Each intersection curve should have multiple points"