Add method to check if point is inside cell (#2905)

* add cell_is_point_inside

* Better parameter description

Co-authored-by: Tetsuo Koyama <tkoyama010@gmail.com>

* Better naming

Co-authored-by: Tetsuo Koyama <tkoyama010@gmail.com>

* 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 <adeak@users.noreply.github.com>

* 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 <adeak@users.noreply.github.com>

* Use typing like syntax in docstring

* black

* raise when vtk returns -1

Co-authored-by: Tetsuo Koyama <tkoyama010@gmail.com>
Co-authored-by: Andras Deak <adeak@users.noreply.github.com>

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

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

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.')

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))

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)