[geom] Implement Mesh.lies_inside(), signed distance

phi/field/_resample.py

@@ -408,19 +408,11 @@ def sample_mesh(f: Field,
                 location: Tensor,
                 gradient: Union[str, Field] = 'green-gauss',
                 order=2,
-                max_steps=2):  # at least 2 to resolve locations outside the mesh
+                max_steps=None):  # at least 2 to resolve locations outside the mesh
+    max_steps = f.mesh._max_cell_walk if max_steps is None else max_steps
     idx = math.find_closest(f.center, location)
     for i in range(max_steps):
-        is_last_step = i == max_steps - 1
-        # --- check if inside, else move to neighbor cell ---
-        closest_normals = f.mesh.face_normals[idx]
-        closest_face_centers = f.mesh.face_centers[idx]
-        offsets = closest_normals.vector @ closest_face_centers.vector  # this dot product could be cashed in the mesh
-        distances = closest_normals.vector @ location.vector - offsets
-        is_outside = math.any(distances > 0, dual)
-        next_idx = math.argmax(distances, dual).index[0]
-        leaves_mesh = next_idx >= instance(f).volume
-        idx = math.where(is_outside & (~leaves_mesh | is_last_step), next_idx, idx)
+        idx, leaves_mesh, is_outside, *_ = f.mesh.cell_walk_towards(location, idx, allow_exit=i == max_steps - 1)
     is_outside_mesh = leaves_mesh & is_outside
     if order <= 1:
         values = rename_dims(f.mesh.pad_boundary(f.values, mode=f.boundary), dual, instance(f))

phi/geom/_mesh.py

@@ -38,7 +38,8 @@ def __init__(self, vertices: Graph,
                  volume: Tensor,
                  faces: Face,
                  valid_mask: Tensor,
-                 face_vertices: Tensor):
+                 face_vertices: Tensor,
+                 max_cell_walk: int = 2):
         """
         Args:
             vertices: Vertex positions, shape (vertices:i, vector:c)
@@ -71,6 +72,7 @@ def __init__(self, vertices: Graph,
         boundary_deltas = (self.face_centers - self.center)[self.all_boundary_faces]
         assert (math.vec_length(boundary_deltas) > 0).all, f"All boundary faces must be separated from the cell centers but 0 distance at the following {channel(math.stored_indices(boundary_deltas)).item_names[0]}:\n{math.nonzero(math.vec_length(boundary_deltas) == 0):full}"
         self._neighbor_offsets = math.concat([cell_deltas, boundary_deltas], '~neighbors')
+        self._max_cell_walk = max_cell_walk
         # --- skewness ---
         # theta_e = math.PI * (vertex_count - 2) / vertex_count
         # e_face =
@@ -225,11 +227,51 @@ def volume(self) -> Tensor:
         return self._volume
 
     def lies_inside(self, location: Tensor) -> Tensor:
-        raise NotImplementedError
+        idx = math.find_closest(self._center, location)
+        for i in range(self._max_cell_walk):
+            idx, leaves_mesh, is_outside, *_ = self.cell_walk_towards(location, idx, allow_exit=i == self._max_cell_walk - 1)
+        return ~(leaves_mesh & is_outside)
 
     def approximate_signed_distance(self, location: Union[Tensor, tuple]) -> Tensor:
+        idx = math.find_closest(self._center, location)
+        for i in range(self._max_cell_walk):
+            idx, leaves_mesh, is_outside, distances, nb_idx = self.cell_walk_towards(location, idx, allow_exit=False)
+        return math.max(distances, dual)
+
+    def approximate_closest_surface(self, location: Tensor) -> Tuple[Tensor, Tensor, Tensor, Tensor, Tensor]:
+        # idx = math.find_closest(self._center, location)
+        # for i in range(self._max_cell_walk):
+        #     idx, leaves_mesh, is_outside, distances, nb_idx = self.cell_walk_towards(location, idx, allow_exit=False)
+        # sgn_dist = math.max(distances, dual)
+        # cell_normals = self.face_normals[idx]
+        # normal = cell_normals[{dual: nb_idx}]
+        # return sgn_dist, delta, normal, offset, face_index
         raise NotImplementedError
 
+    def cell_walk_towards(self, location: Tensor, start_cell_idx: Tensor, allow_exit=False):
+        """
+        If `location` is not within the cell at index `from_cell_idx`, moves to a closer neighbor cell.
+
+        Args:
+            location: Target location as `Tensor`.
+            start_cell_idx: Index of starting cell. Must be a valid cell index.
+            allow_exit: If `True`, returns an invalid index for points outside the mesh, otherwise keeps the current index.
+
+        Returns:
+            index: Index of the neighbor cell or starting cell.
+            leaves_mesh: Whether the walk crossed the mesh boundary. Then `index` is invalid. This is only possible if `allow_exit` is true.
+            is_outside: Whether `location` was outside the cell at index `start_cell_idx`.
+        """
+        closest_normals = self.face_normals[start_cell_idx]
+        closest_face_centers = self.face_centers[start_cell_idx]
+        offsets = closest_normals.vector @ closest_face_centers.vector  # this dot product could be cashed in the mesh
+        distances = closest_normals.vector @ location.vector - offsets
+        is_outside = math.any(distances > 0, dual)
+        nb_idx = math.argmax(distances, dual).index[0]  # cell index or boundary face index
+        leaves_mesh = nb_idx >= instance(self).volume
+        next_idx = math.where(is_outside & (~leaves_mesh | allow_exit), nb_idx, start_cell_idx)
+        return next_idx, leaves_mesh, is_outside, distances, nb_idx
+
     def sample_uniform(self, *shape: math.Shape) -> Tensor:
         raise NotImplementedError
 

tests/commit/geom/test__mesh.py

@@ -1,8 +1,8 @@
 from unittest import TestCase
 
 from phi import math
-from phi.geom import Box, build_mesh, Sphere
-from phiml.math import spatial
+from phi.geom import Box, build_mesh, Sphere, mesh_from_numpy
+from phiml.math import spatial, vec
 
 
 class TestGrid(TestCase):
@@ -14,3 +14,23 @@ def test_build_mesh(self):
         # --- with obstacles ---
         obs = Sphere(x=.5, y=0, radius=.5)
         build_mesh(Box(x=1, y=1), x=10, y=10, obstacles=obs)
+
+    def test_lies_inside(self):
+        points = [(0, 0), (1, 0), (0, 1)]
+        polygons = [(0, 1, 2)]
+        boundaries = {'outer': [(0, 1), (1, 2), (2, 0)]}
+        mesh = mesh_from_numpy(points, polygons, boundaries)
+        math.assert_close(True, mesh.lies_inside(vec(x=.4, y=.4)))
+        math.assert_close(False, mesh.lies_inside(vec(x=1, y=1)))
+        math.assert_close(False, mesh.lies_inside(vec(x=-.1, y=.5)))
+        math.assert_close([False, True], mesh.lies_inside(vec(x=[-.1, .1], y=.5)))
+
+    def test_closest_distance(self):
+        points = [(0, 0), (1, 0), (0, 1)]
+        polygons = [(0, 1, 2)]
+        boundaries = {'outer': [(0, 1), (1, 2), (2, 0)]}
+        mesh = mesh_from_numpy(points, polygons, boundaries)
+        math.assert_close(.1, mesh.approximate_signed_distance(vec(x=-.1, y=.5)))
+        math.assert_close(-.1, mesh.approximate_signed_distance(vec(x=.1, y=.5)))
+        math.assert_close(-.1, mesh.approximate_signed_distance(vec(x=.5, y=.1)))
+        math.assert_close(.1, mesh.approximate_signed_distance(vec(x=.5, y=-.1)))