From d44db9a4d5072578b32904fe3b462bc30fee91bf Mon Sep 17 00:00:00 2001 From: Philipp Holl Date: Mon, 3 Jun 2024 20:50:02 +0200 Subject: [PATCH] [geom] Implement Mesh.lies_inside(), signed distance --- phi/field/_resample.py | 14 +++------- phi/geom/_mesh.py | 46 +++++++++++++++++++++++++++++++-- tests/commit/geom/test__mesh.py | 24 +++++++++++++++-- 3 files changed, 69 insertions(+), 15 deletions(-) diff --git a/phi/field/_resample.py b/phi/field/_resample.py index 35f591393..3820b8b2e 100644 --- a/phi/field/_resample.py +++ b/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)) diff --git a/phi/geom/_mesh.py b/phi/geom/_mesh.py index 81c7295cb..a108d50bb 100644 --- a/phi/geom/_mesh.py +++ b/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 diff --git a/tests/commit/geom/test__mesh.py b/tests/commit/geom/test__mesh.py index d81b66722..c1db5760d 100644 --- a/tests/commit/geom/test__mesh.py +++ b/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)))