Adapt implementation from lace

README.md

@@ -13,6 +13,8 @@ This library is optimized for cloud computation, however it depends on both
 [NumPy][] and [the SciPy k-d tree][ckdtree]. It's designed for use with
 [lacecore][].
 
+Currently works only with triangular meshes.
+
 Prior to version 1.0, this was a library for rendering mesh cross sections to
 SVG. That library has been renamed to [hobart-svg][].
 

hobart/__init__.py

@@ -1 +1,5 @@
+from . import _core
+from ._core import *  # noqa: F403,F401
 from .package_version import __version__  # noqa: F401
+
+__all__ = _core.__all__

hobart/_core.py

@@ -0,0 +1,167 @@
+# Adapted from
+# https://github.com/lace/lace/blob/a19d14eeca37a66af442c5a73c3934987d2d7d9a/lace/intersection.py
+# https://github.com/lace/lace/blob/a19d14eeca37a66af442c5a73c3934987d2d7d9a/lace/test_intersection.py
+# By Body Labs & Metabolize, BSD License
+
+import numpy as np
+from polliwog import Polyline
+from scipy.spatial import cKDTree
+import vg
+from ._internal import EdgeMap, Graph
+from ._validation import check_indices
+
+__all__ = ["faces_intersecting_plane", "intersect_mesh_with_plane"]
+
+
+def faces_intersecting_plane(vertices, faces, plane):
+    """
+    Efficiently compute which of the given faces intersect the given plane.
+
+    Args:
+        vertices (np.ndarray): The vertices, as a `nx3` array.
+        faces (np.ndarray): The face indices, as a `kx3` array.
+        plane (polliwog.Plane): The plane of interest.
+
+    Returns:
+        np.ndarray: A boolean mask indicating the faces which intersect
+            the given plane.
+
+    See also:
+        https://polliwog.readthedocs.io/en/latest/#polliwog.Plane
+    """
+    vg.shape.check(locals(), "vertices", (-1, 3))
+    vg.shape.check(locals(), "faces", (-1, 3))
+    check_indices(faces, len(vertices), "faces")
+
+    signed_distances = plane.signed_distance(vertices)
+    return np.abs(np.sign(signed_distances)[faces].sum(axis=1)) != 3
+
+
+def intersect_mesh_with_plane(
+    vertices, faces, plane, neighborhood=None, ret_pointcloud=False
+):
+    """
+    Takes a cross section of planar point cloud with a Mesh object.
+    Ignore those points which intersect at a vertex - the probability of
+    this event is small, and accounting for it complicates the algorithm.
+
+    When a plane may intersect the mesh more than once, provide a
+    `neighborhood`, which is a set of points. The cross section closest
+    to this list of points is returned (using a kd-tree).
+
+    Args:
+        vertices (np.ndarray): The vertices, as a `nx3` array.
+        faces (np.ndarray): The face indices, as a `kx3` array.
+        plane (polliwog.Plane): The plane of interest.
+        neighborhood (np.ndarray): One or more points of interest, used to
+            select the desired cross section when a plane may intersect more
+            than once.
+        ret_pointcloud (bool): When `True`, return an unstructured point
+            cloud instead of a list of polylines. This is useful when
+            you aren't specifying a neighborhood and you only care about
+            e.g. some apex of the intersection.
+
+    Returns:
+        list: A list of `polliwog.Polyline` instances.
+
+    See also:
+        https://polliwog.readthedocs.io/en/latest/#polliwog.Plane
+        https://polliwog.readthedocs.io/en/latest/#polliwog.Polyline
+    """
+    if neighborhood is not None:
+        vg.shape.check(locals(), "neighborhood", (-1, 3))
+
+    # 1: Select those faces that intersect the plane, fs
+    fs = faces[faces_intersecting_plane(vertices, faces, plane)]
+
+    if len(fs) == 0:
+        # Nothing intersects.
+        if ret_pointcloud:
+            return np.zeros((0, 3))
+        elif neighborhood is not None:
+            return None
+        else:
+            return []
+
+    # and edges of those faces
+    es = np.vstack((fs[:, (0, 1)], fs[:, (1, 2)], fs[:, (2, 0)]))
+
+    # 2: Find the edges where each of those faces actually cross the plane
+    intersection_map = EdgeMap()
+
+    pts, pt_is_valid = plane.line_segment_xsections(
+        vertices[es[:, 0]], vertices[es[:, 1]]
+    )
+    valid_pts = pts[pt_is_valid]
+    valid_es = es[pt_is_valid]
+    for val, e in zip(valid_pts, valid_es):
+        if not intersection_map.contains(e[0], e[1]):
+            intersection_map.add(e[0], e[1], val)
+    verts = np.array(intersection_map.values)
+
+    # 3: Build the edge adjacency graph
+    G = Graph(verts.shape[0])
+    for f in fs:
+        # Since we're dealing with a triangle that intersects the plane,
+        # exactly two of the edges will intersect (note that the only other
+        # sorts of "intersections" are one edge in plane or all three edges in
+        # plane, which won't be picked up by mesh_intersecting_faces).
+        e0 = intersection_map.index(f[0], f[1])
+        e1 = intersection_map.index(f[0], f[2])
+        e2 = intersection_map.index(f[1], f[2])
+        if e0 is None:
+            G.add_edge(e1, e2)
+        elif e1 is None:
+            G.add_edge(e0, e2)
+        else:
+            G.add_edge(e0, e1)  # pragma: no cover
+
+    # 4: Find the paths for each component
+    components = []
+    components_closed = []
+    while len(G) > 0:
+        path = G.pop_euler_path()
+        if path is None:
+            raise ValueError(  # pragma: no cover
+                "Mesh slice has too many odd degree edges; can't find a path along the edge"
+            )
+        component_verts = verts[path]
+
+        if np.all(component_verts[0] == component_verts[-1]):
+            # Because the closed polyline will make that last link:
+            component_verts = np.delete(component_verts, 0, axis=0)
+            components_closed.append(True)
+        else:
+            components_closed.append(False)
+        components.append(component_verts)
+
+    # 6 (optional - only if 'neighborhood' is provided): Use a KDTree to select
+    # the component with minimal distance to 'neighborhood'.
+    if neighborhood is not None and len(components) > 1:
+        tree = cKDTree(neighborhood)
+
+        # The number of components will not be large in practice, so this loop
+        # won't hurt.
+        means = [np.mean(tree.query(component)[0]) for component in components]
+        index = np.argmin(means)
+        if ret_pointcloud:
+            return components[index]
+        else:
+            return Polyline(components[index], is_closed=components_closed[index])
+    elif neighborhood is not None and len(components) == 1:
+        if ret_pointcloud:  # pragma: no cover
+            return components[0]  # pragma: no cover
+        else:
+            return Polyline(
+                components[0], is_closed=components_closed[0]
+            )  # pragma: no cover
+    else:
+        # No neighborhood provided, so return all the components, either in a
+        # pointcloud or as separate polylines.
+        if ret_pointcloud:
+            return np.vstack(components)
+        else:
+            return [
+                Polyline(v, is_closed=closed)
+                for v, closed in zip(components, components_closed)
+            ]

hobart/_internal.py

@@ -0,0 +1,95 @@
+class EdgeMap:
+    """
+    A quick two-level dictionary where the two keys are interchangeable (i.e.
+    a symmetric graph).
+    """
+
+    def __init__(self):
+        # Store indicies into self.values here, to make it easier to get inds
+        # or values.
+        self.d = {}
+        self.values = []
+
+    def _order(self, u, v):
+        if u < v:
+            return u, v
+        else:
+            return v, u
+
+    def add(self, u, v, val):
+        low, high = self._order(u, v)
+        if low not in self.d:
+            self.d[low] = {}
+        self.values.append(val)
+        self.d[low][high] = len(self.values) - 1
+
+    def contains(self, u, v):
+        low, high = self._order(u, v)
+        if low in self.d and high in self.d[low]:
+            return True
+        return False
+
+    def index(self, u, v):
+        low, high = self._order(u, v)
+        try:
+            return self.d[low][high]
+        except KeyError:
+            return None
+
+
+class Graph:
+    """
+    A little utility class to build a symmetric graph and calculate Euler Paths.
+    """
+
+    def __init__(self, size):
+        self.size = size
+        self.d = {}
+
+    def __len__(self):
+        return len(self.d)
+
+    def add_edge(self, u, v):
+        assert u >= 0 and u < self.size
+        assert v >= 0 and v < self.size
+        if u not in self.d:
+            self.d[u] = set()
+        if v not in self.d:
+            self.d[v] = set()
+        self.d[u].add(v)
+        self.d[v].add(u)
+
+    def remove_edge(self, u, v):
+        if u in self.d and v in self.d[u]:
+            self.d[u].remove(v)
+        if v in self.d and u in self.d[v]:
+            self.d[v].remove(u)  # pragma: no cover
+        if v in self.d and len(self.d[v]) == 0:
+            del self.d[v]
+        if u in self.d and len(self.d[u]) == 0:
+            del self.d[u]
+
+    def pop_euler_path(self, allow_multiple_connected_components=True):
+        # Based on code from Przemek Drochomirecki, Krakow, 5 Nov 2006
+        # http://code.activestate.com/recipes/498243-finding-eulerian-path-in-undirected-graph/
+        # Under PSF License
+        # NB: MUTATES d
+
+        # Count the number of vertices with odd degree.
+        odd = [x for x in self.d if len(self.d[x]) & 1]
+        odd.append(list(self.d.keys())[0])
+        if not allow_multiple_connected_components and len(odd) > 3:
+            return None  # pragma: no cover
+        stack = [odd[0]]
+        path = []
+
+        # Main algorithm.
+        while stack:
+            v = stack[-1]
+            if v in self.d:
+                u = self.d[v].pop()
+                stack.append(u)
+                self.remove_edge(u, v)
+            else:
+                path.append(stack.pop())
+        return path

hobart/_validation.py

@@ -0,0 +1,9 @@
+import numpy as np
+
+
+def check_indices(indices, num_elements, name):
+    # Vendored in from lacecore.
+    if np.any(indices >= num_elements):
+        raise ValueError(
+            "Expected indices in {} to be less than {}".format(name, num_elements)
+        )