Fix point on plane detection for NumPy on Linux

Closes #198

polliwog/plane/_plane_functions.py

@@ -4,6 +4,7 @@
 from .._common.shape import check_shape_any, columnize
 
 __all__ = [
+    "EPSILON_COPLANAR",
     "plane_normal_from_points",
     "plane_equation_from_points",
     "normal_and_offset_from_plane_equations",
@@ -12,6 +13,8 @@
     "mirror_point_across_plane",
 ]
 
+EPSILON_COPLANAR = 1e-12
+
 
 def plane_normal_from_points(points, normalize=True):
     """
@@ -62,20 +65,38 @@ def normal_and_offset_from_plane_equations(plane_equations):
     return normal, offset
 
 
-def signed_distance_to_plane(points, plane_equations):
+def signed_distance_to_plane(points, plane_equations, epsilon=EPSILON_COPLANAR):
     """
     Return the signed distances from each point to the corresponding plane.
 
     For convenience, can also be called with a single point and a single
     plane.
+
+    Args:
+        points (np.ndarray): The points of interest as `kx3`.
+        plane_equations (np.ndarray): The plane equations as `kx4`.
+        epsilon (float): Return 0 for points within this distance of the
+            plane. Pass `None` to skip this check.
+
+    Return:
+        The signed distance, or for stacked inputs, an array of signed
+            distances.
     """
     k = check_shape_any(points, (3,), (-1, 3), name="points")
     check_shape_any(
         plane_equations, (4,), (-1 if k is None else k, 4), name="plane_equations"
     )
 
     normals, offsets = normal_and_offset_from_plane_equations(plane_equations)
-    return vg.dot(points, normals) + offsets
+    result = vg.dot(points, normals) + offsets
+
+    if epsilon is not None:
+        if np.isscalar(result):
+            if np.abs(result) < epsilon:
+                result = 0.0
+        else:
+            result[np.abs(result) < epsilon] = 0.0
+    return result
 
 
 def translate_points_along_plane_normal(points, plane_equations, factor):
@@ -95,7 +116,7 @@ def translate_points_along_plane_normal(points, plane_equations, factor):
     assert isinstance(factor, numbers.Real)
 
     # Translate the point back to the plane along the normal.
-    signed_distance = signed_distance_to_plane(points, plane_equations)
+    signed_distance = signed_distance_to_plane(points, plane_equations, epsilon=None)
     normals, _ = normal_and_offset_from_plane_equations(plane_equations)
 
     if np.isscalar(signed_distance):

polliwog/plane/_plane_object.py

@@ -1,6 +1,7 @@
 import numpy as np
 import vg
 from ._plane_functions import (
+    EPSILON_COPLANAR,
     mirror_point_across_plane,
     plane_normal_from_points,
     project_point_to_plane,
@@ -138,14 +139,25 @@ def flipped(self):
         """
         return Plane(point_on_plane=self._r0, unit_normal=-self._n)
 
-    def sign(self, points):
+    def sign(self, points, epsilon=EPSILON_COPLANAR):
         """
         Given an array of points, return an array with +1 for points in front
         of the plane (in the direction of the normal), -1 for points behind
         the plane (away from the normal), and 0 for points on the plane.
 
+        Args:
+            points (np.arraylike): A 3D point or a `kx3` stack of points.
+            epsilon (float): Return 0 for points within this distance of the
+                plane.
+
+        Returns:
+            depends:
+
+            - Given a single 3D point, the distance as a NumPy scalar.
+            - Given a `kx3` stack of points, an `k` array of distances.
+
         """
-        return np.sign(self.signed_distance(points))
+        return np.sign(self.signed_distance(points, epsilon=epsilon))
 
     def points_in_front(self, points, inverted=False, ret_indices=False):
         """
@@ -203,21 +215,23 @@ def points_on_or_in_front(self, points, inverted=False, ret_indices=False):
 
         return indices if ret_indices else points[indices]
 
-    def signed_distance(self, points):
+    def signed_distance(self, points, epsilon=EPSILON_COPLANAR):
         """
         Returns the signed distances to the given points or the signed
         distance to a single point.
 
         Args:
             points (np.arraylike): A 3D point or a `kx3` stack of points.
+            epsilon (float): Return 0 for points within this distance of the
+                plane.
 
         Returns:
             depends:
 
             - Given a single 3D point, the distance as a NumPy scalar.
             - Given a `kx3` stack of points, an `k` array of distances.
         """
-        return signed_distance_to_plane(points, self.equation)
+        return signed_distance_to_plane(points, self.equation, epsilon=epsilon)
 
     def distance(self, points):
         return np.absolute(self.signed_distance(points))

polliwog/polyline/_polyline_object.py

@@ -411,7 +411,7 @@ def intersect_plane(self, plane, ret_edge_indices=False):
         """
         # TODO: Refactor to use `..plane.intersections.intersect_segment_with_plane()`.
         # Identify edges with endpoints that are not on the same side of the plane
-        signed_distances = plane.signed_distance(self.v)
+        signed_distances = plane.signed_distance(self.v, epsilon=None)
         which_es = np.abs(np.sign(signed_distances)[self.e].sum(axis=1)) != 2
         # For the intersecting edges, compute the distance of the endpoints to the plane
         endpoint_distances = np.abs(signed_distances[self.e[which_es]])

requirements.txt

@@ -1,3 +1,3 @@
-numpy<1.19.0
+numpy
 ounce>=1.1.0,<2.0
 vg>=1.5.0

requirements_doc.txt

@@ -1,5 +1,5 @@
 # Imported from code.
-numpy<1.19.0
+numpy
 ounce>=1.1.0,<2.0
 vg>=1.5.0