From fdc66f3735717a3fc72275044080577112861921 Mon Sep 17 00:00:00 2001 From: MatthewFlamm <39341281+MatthewFlamm@users.noreply.github.com> Date: Thu, 7 Jul 2022 17:56:25 -0400 Subject: [PATCH] Add method to check if point is inside cell (#2905) * add cell_is_point_inside * Better parameter description Co-authored-by: Tetsuo Koyama * Better naming Co-authored-by: Tetsuo Koyama * use better name in doctest and test * fix for python 3.10 * use collections.abc.Sequence for all Python versions * Use interval check Co-authored-by: Andras Deak * allow points * revert unintended removal of whitespace * check for int type for ind * return 2d array when a 2d array is supplied * return ndarray instead of list * whitespace * Fix unintended emacs closing file problem * add cell_point_ids to avoid painful diff merge * fix whitespace in _coerce_pointslike_arg * fix versionchanged format * Some minor docstring wording/markup changes * Remove return and and very description for second return * Also use numpy.integer Co-authored-by: Andras Deak * Use typing like syntax in docstring * black * raise when vtk returns -1 Co-authored-by: Tetsuo Koyama Co-authored-by: Andras Deak --- pyvista/core/dataset.py | 88 +++++++++++++++++++++++--- pyvista/utilities/common.py | 14 ++-- pyvista/utilities/geometric_objects.py | 2 +- tests/test_dataset.py | 28 ++++++++ tests/utilities/test_common.py | 31 ++++++--- 5 files changed, 140 insertions(+), 23 deletions(-) diff --git a/pyvista/core/dataset.py b/pyvista/core/dataset.py index 43bd5c9b2f..b3f769e3ee 100644 --- a/pyvista/core/dataset.py +++ b/pyvista/core/dataset.py @@ -416,7 +416,7 @@ def points(self, points: Union[VectorArray, NumericArray, _vtk.vtkPoints]): self.Modified() return # otherwise, wrap and use the array - points = _coerce_pointslike_arg(points, copy=False) + points, _ = _coerce_pointslike_arg(points, copy=False) vtk_points = pyvista.vtk_points(points, False) if not pdata: self.SetPoints(vtk_points) @@ -2187,11 +2187,19 @@ def find_closest_cell( Index or indices of the cell in this mesh that is/are closest to the given point(s). + .. versionchanged:: 0.35.0 + Inputs of shape ``(1, 3)`` now return a :class:`numpy.ndarray` + of shape ``(1,)``. + numpy.ndarray Point or points inside a cell of the mesh that is/are closest to the given point(s). Only returned if ``return_closest_point=True``. + .. versionchanged:: 0.35.0 + Inputs of shape ``(1, 3)`` now return a :class:`numpy.ndarray` + of the same shape. + Warnings -------- This method may still return a valid cell index even if the point @@ -2260,7 +2268,7 @@ def find_closest_cell( array([1., 1., 0.]) """ - point = _coerce_pointslike_arg(point, copy=False) + point, singular = _coerce_pointslike_arg(point, copy=False) locator = _vtk.vtkCellLocator() locator.SetDataSet(self) @@ -2282,11 +2290,9 @@ def find_closest_cell( closest_points.append(closest_point) out_cells: Union[int, np.ndarray] = ( - closest_cells[0] if len(closest_cells) == 1 else np.array(closest_cells) - ) - out_points = ( - np.array(closest_points[0]) if len(closest_points) == 1 else np.array(closest_points) + closest_cells[0] if singular else np.array(closest_cells) ) + out_points = np.array(closest_points[0]) if singular else np.array(closest_points) if return_closest_point: return out_cells, out_points @@ -2309,6 +2315,10 @@ def find_containing_cell( Index or indices of the cell in this mesh that contains the given point. + .. versionchanged:: 0.35.0 + Inputs of shape ``(1, 3)`` now return a :class:`numpy.ndarray` + of shape ``(1,)``. + See Also -------- DataSet.find_closest_point @@ -2342,14 +2352,14 @@ def find_containing_cell( (1000,) """ - point = _coerce_pointslike_arg(point, copy=False) + point, singular = _coerce_pointslike_arg(point, copy=False) locator = _vtk.vtkCellLocator() locator.SetDataSet(self) locator.BuildLocator() containing_cells = [locator.FindCell(node) for node in point] - return containing_cells[0] if len(containing_cells) == 1 else np.array(containing_cells) + return containing_cells[0] if singular else np.array(containing_cells) def find_cells_along_line( self, @@ -2568,3 +2578,65 @@ def cell_point_ids(self, ind: int) -> List[int]: cell = self.GetCell(ind) point_ids = cell.GetPointIds() return [point_ids.GetId(i) for i in range(point_ids.GetNumberOfIds())] + + def point_is_inside_cell( + self, ind: int, point: Union[VectorArray, NumericArray] + ) -> Union[int, np.ndarray]: + """Return whether one or more points are inside a cell. + + .. versionadded:: 0.35.0 + + Parameters + ---------- + ind : int + Cell ID. + + point : Sequence[float] or np.ndarray + Coordinates of point to query (length 3) or a ``numpy`` array of ``n`` + points with shape ``(n, 3)``. + + Returns + ------- + bool or numpy.ndarray + Whether point(s) is/are inside cell. A scalar bool is only returned if + the input point has shape ``(3,)``. + + Examples + -------- + >>> from pyvista import examples + >>> mesh = examples.load_hexbeam() + >>> mesh.cell_bounds(0) + (0.0, 0.5, 0.0, 0.5, 0.0, 0.5) + >>> mesh.point_is_inside_cell(0, [0.2, 0.2, 0.2]) + True + + """ + if not isinstance(ind, (int, np.integer)): + raise TypeError(f"ind must be an int, got {type(ind)}") + + if not 0 <= ind < self.n_cells: + raise ValueError(f"ind must be >= 0 and < {self.n_cells}, got {ind}") + + co_point, singular = _coerce_pointslike_arg(point, copy=False) + + cell = self.GetCell(ind) + npoints = cell.GetPoints().GetNumberOfPoints() + + closest_point = [0.0, 0.0, 0.0] + sub_id = _vtk.mutable(0) + pcoords = [0.0, 0.0, 0.0] + dist2 = _vtk.mutable(0.0) + weights = [0.0] * npoints + + in_cell = np.empty(shape=co_point.shape[0], dtype=np.bool_) + for i, node in enumerate(co_point): + is_inside = cell.EvaluatePosition(node, closest_point, sub_id, pcoords, dist2, weights) + if not 0 <= is_inside <= 1: + raise RuntimeError( + f"Computational difficulty encountered for point {node} in cell {ind}" + ) + in_cell[i] = bool(is_inside) + + if singular: + return in_cell[0] + return in_cell diff --git a/pyvista/utilities/common.py b/pyvista/utilities/common.py index 5f894d51e8..40f92df8b4 100644 --- a/pyvista/utilities/common.py +++ b/pyvista/utilities/common.py @@ -1,6 +1,6 @@ """Common functions.""" import collections.abc -from typing import Union +from typing import Tuple, Union import numpy as np @@ -9,7 +9,7 @@ def _coerce_pointslike_arg( points: Union[NumericArray, VectorArray], copy: bool = False -) -> np.ndarray: +) -> Tuple[np.ndarray, bool]: """Check and coerce arg to (n, 3) np.ndarray. Parameters @@ -25,6 +25,8 @@ def _coerce_pointslike_arg( ------- np.ndarray Size (n, 3) array. + bool + Whether the input was a single point in an array-like with shape (3,). """ if isinstance(points, collections.abc.Sequence): @@ -35,14 +37,18 @@ def _coerce_pointslike_arg( if points.ndim > 2: raise ValueError("Array of points must be 1D or 2D") + if points.ndim == 2: if points.shape[1] != 3: raise ValueError("Array of points must have three values per point (shape (n, 3))") + singular = False + else: if points.size != 3: raise ValueError("Given point must have three values") + singular = True points = np.reshape(points, [1, 3]) if copy: - return points.copy() - return points + return points.copy(), singular + return points, singular diff --git a/pyvista/utilities/geometric_objects.py b/pyvista/utilities/geometric_objects.py index 0d47403f69..1dc885667c 100644 --- a/pyvista/utilities/geometric_objects.py +++ b/pyvista/utilities/geometric_objects.py @@ -487,7 +487,7 @@ def MultipleLines(points=[[-0.5, 0.0, 0.0], [0.5, 0.0, 0.0]]): >>> mesh = pyvista.MultipleLines(points=[[0, 0, 0], [1, 1, 1], [0, 0, 1]]) >>> mesh.plot(color='k', line_width=10) """ - points = _coerce_pointslike_arg(points) + points, _ = _coerce_pointslike_arg(points) src = _vtk.vtkLineSource() if not (len(points) >= 2): raise ValueError('>=2 points need to define multiple lines.') diff --git a/tests/test_dataset.py b/tests/test_dataset.py index 3ce92d68c8..dea4417750 100644 --- a/tests/test_dataset.py +++ b/tests/test_dataset.py @@ -1102,6 +1102,34 @@ def test_cell_type(grid): assert isinstance(ctype, int) +def test_point_is_inside_cell(): + grid = pyvista.UniformGrid(dims=(2, 2, 2)) + assert grid.point_is_inside_cell(0, [0.5, 0.5, 0.5]) + assert not grid.point_is_inside_cell(0, [-0.5, -0.5, -0.5]) + + assert grid.point_is_inside_cell(0, np.array([0.5, 0.5, 0.5])) + + # cell ind out of range + with pytest.raises(ValueError): + grid.point_is_inside_cell(100000, [0.5, 0.5, 0.5]) + with pytest.raises(ValueError): + grid.point_is_inside_cell(-1, [0.5, 0.5, 0.5]) + + # cell ind wrong type + with pytest.raises(TypeError): + grid.point_is_inside_cell(0.1, [0.5, 0.5, 0.5]) + + # point not well formed + with pytest.raises(TypeError): + grid.point_is_inside_cell(0, 0.5) + with pytest.raises(ValueError): + grid.point_is_inside_cell(0, [0.5, 0.5]) + + # multi-dimensional + in_cell = grid.point_is_inside_cell(0, [[0.5, 0.5, 0.5], [-0.5, -0.5, -0.5]]) + assert np.array_equal(in_cell, np.array([True, False])) + + def test_serialize_deserialize(datasets): for dataset in datasets: dataset_2 = pickle.loads(pickle.dumps(dataset)) diff --git a/tests/utilities/test_common.py b/tests/utilities/test_common.py index a54401395d..86a4e15823 100644 --- a/tests/utilities/test_common.py +++ b/tests/utilities/test_common.py @@ -9,57 +9,68 @@ def test_coerce_point_like_arg(): # Test with Sequence point = [1.0, 2.0, 3.0] - coerced_arg = _coerce_pointslike_arg(point) + coerced_arg, singular = _coerce_pointslike_arg(point) assert isinstance(coerced_arg, np.ndarray) assert coerced_arg.shape == (1, 3) assert np.array_equal(coerced_arg, np.array([[1.0, 2.0, 3.0]])) + assert singular # Test with 1D np.ndarray point = np.array([1.0, 2.0, 3.0]) - coerced_arg = _coerce_pointslike_arg(point) + coerced_arg, singular = _coerce_pointslike_arg(point) assert isinstance(coerced_arg, np.ndarray) assert coerced_arg.shape == (1, 3) assert np.array_equal(coerced_arg, np.array([[1.0, 2.0, 3.0]])) + assert singular + + # Test with (1, 3) np.ndarray + point = np.array([[1.0, 2.0, 3.0]]) + coerced_arg, singular = _coerce_pointslike_arg(point) + assert isinstance(coerced_arg, np.ndarray) + assert coerced_arg.shape == (1, 3) + assert np.array_equal(coerced_arg, np.array([[1.0, 2.0, 3.0]])) + assert not singular # Test with 2D ndarray point = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]) - coerced_arg = _coerce_pointslike_arg(point) + coerced_arg, singular = _coerce_pointslike_arg(point) assert isinstance(coerced_arg, np.ndarray) assert coerced_arg.shape == (2, 3) assert np.array_equal(coerced_arg, np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])) + assert not singular def test_coerce_point_like_arg_copy(): # Sequence is always copied point = [1.0, 2.0, 3.0] - coerced_arg = _coerce_pointslike_arg(point, copy=True) + coerced_arg, _ = _coerce_pointslike_arg(point, copy=True) point[0] = 10.0 assert np.array_equal(coerced_arg, np.array([[1.0, 2.0, 3.0]])) point = [1.0, 2.0, 3.0] - coerced_arg = _coerce_pointslike_arg(point, copy=False) + coerced_arg, _ = _coerce_pointslike_arg(point, copy=False) point[0] = 10.0 assert np.array_equal(coerced_arg, np.array([[1.0, 2.0, 3.0]])) # 1D np.ndarray can be copied or not point = np.array([1.0, 2.0, 3.0]) - coerced_arg = _coerce_pointslike_arg(point, copy=True) + coerced_arg, _ = _coerce_pointslike_arg(point, copy=True) point[0] = 10.0 assert np.array_equal(coerced_arg, np.array([[1.0, 2.0, 3.0]])) point = np.array([1.0, 2.0, 3.0]) - coerced_arg = _coerce_pointslike_arg(point, copy=False) + coerced_arg, _ = _coerce_pointslike_arg(point, copy=False) point[0] = 10.0 assert np.array_equal(coerced_arg, np.array([[10.0, 2.0, 3.0]])) # 2D np.ndarray can be copied or not point = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]) - coerced_arg = _coerce_pointslike_arg(point, copy=True) + coerced_arg, _ = _coerce_pointslike_arg(point, copy=True) point[0, 0] = 10.0 assert np.array_equal(coerced_arg, np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])) point = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]) - coerced_arg = _coerce_pointslike_arg(point, copy=False) + coerced_arg, _ = _coerce_pointslike_arg(point, copy=False) point[0, 0] = 10.0 assert np.array_equal(coerced_arg, np.array([[10.0, 2.0, 3.0], [4.0, 5.0, 6.0]])) @@ -87,7 +98,7 @@ def test_coerce_point_like_arg_errors(): def test_coerce_points_like_args_does_not_copy(): source = np.random.rand(100, 3) - output = _coerce_pointslike_arg(source) # test that copy=False is default + output, _ = _coerce_pointslike_arg(source) # test that copy=False is default output /= 2 assert np.array_equal(output, source) assert np.may_share_memory(output, source)