From 4f433a7e5533ef2fb8355a5faf105fb316fc852c Mon Sep 17 00:00:00 2001 From: Paul Melnikow Date: Fri, 3 Jul 2020 16:53:56 -0400 Subject: [PATCH] Fix point on plane detection for NumPy on Linux Closes #198 --- polliwog/plane/_plane_functions.py | 27 ++++++++++++++++++++++++--- polliwog/plane/_plane_object.py | 22 ++++++++++++++++++---- polliwog/polyline/_polyline_object.py | 2 +- requirements.txt | 2 +- requirements_doc.txt | 2 +- 5 files changed, 45 insertions(+), 10 deletions(-) diff --git a/polliwog/plane/_plane_functions.py b/polliwog/plane/_plane_functions.py index fd2bb10..4de26a0 100644 --- a/polliwog/plane/_plane_functions.py +++ b/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,12 +65,22 @@ 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( @@ -75,7 +88,15 @@ def signed_distance_to_plane(points, 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): diff --git a/polliwog/plane/_plane_object.py b/polliwog/plane/_plane_object.py index 78f07a4..f36d475 100644 --- a/polliwog/plane/_plane_object.py +++ b/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,13 +215,15 @@ 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: @@ -217,7 +231,7 @@ def signed_distance(self, points): - 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)) diff --git a/polliwog/polyline/_polyline_object.py b/polliwog/polyline/_polyline_object.py index c04e4d3..4d35c2f 100644 --- a/polliwog/polyline/_polyline_object.py +++ b/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]]) diff --git a/requirements.txt b/requirements.txt index 0399f29..40dc8ea 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,3 @@ -numpy<1.19.0 +numpy ounce>=1.1.0,<2.0 vg>=1.5.0 diff --git a/requirements_doc.txt b/requirements_doc.txt index 56da9ef..f3406bd 100644 --- a/requirements_doc.txt +++ b/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