Commit 9469dd0e authored by dboe's avatar dboe
Browse files

added test for _in_triangles and in_triangles

parent 0135c749
import numpy as np
import unittest
import tfields
class Triangles3D_Test(unittest.TestCase):
def setUp(self):
self.mesh = tfields.Mesh3D([[1, 0, 0], [0, 1, 0], [0, 0, 0]], faces=[[0, 1, 2]])
self.inst = self.mesh.triangles()
self.point = tfields.Points3D([[0.2, 0.2, 0]])
self.expected_mask = np.array([True], dtype=bool)
self.points = tfields.Points3D([[0.2, 0.2, 0], [0.2, 0.2, 0]])
self.expected_face_indices = np.array([0, 0], dtype=int)
self.delta = 0.0
def test__in_triangles(self):
self.assertTrue(
np.array_equal(
self.inst._in_triangles(self.point, delta=self.delta),
self.expected_mask,
)
)
def test_in_triangles(self):
indices = self.inst.in_triangles(self.points, delta=self.delta)
self.assertTrue(np.array_equal(indices, self.expected_face_indices))
class Triangles3D_Test_Raises(Triangles3D_Test):
def test_raises(self):
# Raises:
# Wrong format of point will raise a ValueError
with self.assertRaises(ValueError) as ctx:
self.inst._in_triangles(self.points)
self.assertEqual(
"point must be castable to shape 3 but is of shape (2, 3)",
str(ctx.exception),
)
def test__in_triangles_not_invertible(self):
# If you define triangles that have colinear side vectors or in general lead to
# not invertable matrices the you will always get False
mesh = tfields.Mesh3D(
[[0, 0, 0], [2, 0, 0], [4, 0, 0], [0, 1, 0]], faces=[[0, 1, 2], [0, 1, 3]]
)
with self.assertWarns(UserWarning):
mask = mesh.triangles()._in_triangles(np.array([0.2, 0.2, 0]), delta=0.3)
self.assertTrue(np.array_equal(mask, np.array([False, True], dtype=bool)))
class Triangles3D_Test_TwoFaces(Triangles3D_Test):
def setUp(self):
self.mesh = tfields.Mesh3D(
[[1, 0, 0], [0, 1, 0], [0, 0, 0], [4, 0, 0], [4, 4, 0], [8, 0, 0]],
faces=[[0, 1, 2], [3, 4, 5]],
)
self.inst = self.mesh.triangles()
self.point = np.array([0.2, 0.2, 0])
self.expected_mask = np.array([True, False], dtype=bool)
self.points = tfields.Points3D([self.point])
self.expected_face_indices = np.array([0], dtype=int)
self.delta = 0.0
class Triangles3D_Test_TwoFacesOther(Triangles3D_Test_TwoFaces):
def setUp(self):
super().setUp()
self.point = np.array([5, 2, 0])
self.expected_mask = np.array([False, True], dtype=bool)
self.points = tfields.Points3D([self.point])
self.expected_face_indices = np.array([1], dtype=int)
class Triangles3D_TestDeltaOut(Triangles3D_Test):
def setUp(self):
self.mesh = tfields.Mesh3D(
[[1, 0, 0], [0, 1, 0], [0, 0, 0], [4, 0, 0], [4, 4, 0], [8, 0, 0]],
faces=[[0, 1, 2], [3, 4, 5]],
)
self.inst = self.mesh.triangles()
self.point = np.array([0.2, 0.2, 9000])
self.expected_mask = np.array([False, False], dtype=bool)
self.points = tfields.Points3D([self.point])
self.expected_face_indices = np.array([-1], dtype=int)
self.delta = 0.0
class Triangles3D_TestDeltaIn(Triangles3D_TestDeltaOut):
def setUp(self):
super().setUp()
self.point = np.array([0.2, 0.2, 0.1])
self.expected_mask = np.array([True, False], dtype=bool)
self.points = tfields.Points3D([self.point])
self.expected_face_indices = np.array([0], dtype=int)
self.delta = 0.2
class Triangles3D_TestDeltaNone(Triangles3D_TestDeltaOut):
# if you set delta to None, the minimal distance point(s) are accepted automatically
def setUp(self):
super().setUp()
self.point = np.array([0.2, 0.2, 0.1])
self.expected_mask = np.array([True, False], dtype=bool)
self.points = tfields.Points3D([self.point])
self.expected_face_indices = np.array([0], dtype=int)
self.delta = None
......@@ -451,20 +451,16 @@ class Searcher(Node):
>>> assert np.array_equal(box_res, usual_res)
"""
# raise ValueError("Broken feature. We are working on it!")
if not self.is_root():
raise ValueError("in_faces may only be called by root Node.")
if self.at_intersection != "keep":
raise ValueError(
"Intersection method must be 'keep' for in_faces" "method."
)
if self.mesh.nfaces() == 0:
return np.empty((tensors.shape[0], 0), dtype=bool)
assert self.at_intersection == "keep", "Intersection method must be 'keep'"
if delta == -1:
delta = self.delta
masks = np.zeros((len(tensors), self.mesh.nfaces()), dtype=bool)
indices = np.full(tensors.shape[0], np.nan, dtype=int)
if self.mesh.nfaces() == 0:
return indices
with tensors.tmp_transform(self.mesh.coord_sys):
for i, point in enumerate(iter(tensors)):
leaf = self.find_leaf(point)
......@@ -476,8 +472,8 @@ class Searcher(Node):
original_face_indices = leaf.template.maps[3].fields[0][leaf_mask]
if not assign_multiple and len(original_face_indices) > 0:
original_face_indices = original_face_indices[:1]
masks[i, original_face_indices] = True
return masks
indices[i] = original_face_indices
return indices
if __name__ == "__main__":
......
......@@ -533,25 +533,21 @@ class Mesh3D(tfields.TensorMaps):
def nfaces(self):
return self.faces.shape[0]
def in_faces(self, points, delta, assign_multiple=False):
def in_faces(self, points, delta, **kwargs):
"""
Check whether points lie within triangles with Barycentric Technique
see Triangles3D.in_triangles
If multiple requests are done on huge meshes,
see Triangles3D.in_triangles. If multiple requests are done on huge meshes,
this can be hugely optimized by requesting the property
self.tree or setting it to self.tree = <saved tree> before
calling in_faces
self.tree or setting it to self.tree = <saved tree> before calling in_faces.
"""
key = "mesh_tree"
if hasattr(self, "_cache") and key in self._cache:
log = logging.getLogger()
log.info("Using cached decision tree to speed up point - face mapping.")
masks = self.tree.in_faces(points, delta, assign_multiple=assign_multiple)
indices = self.tree.in_faces(points, delta, **kwargs)
else:
masks = self.triangles().in_triangles(
points, delta, assign_multiple=assign_multiple
)
return masks
indices = self.triangles().in_triangles(points, delta, **kwargs)
return indices
@property
def tree(self):
......@@ -560,15 +556,14 @@ class Mesh3D(tfields.TensorMaps):
hugely optimize 'in_faces' searches
Examples:
# >>> mesh = tfields.Mesh3D.grid((0, 1, 3), (1, 2, 3), (2, 3, 3))
# >>> _ = mesh.tree
# >>> assert hasattr(mesh, '_cache')
# >>> assert 'mesh_tree' in mesh._cache
# >>> mask = mesh.in_faces(tfields.Points3D([[0.2, 1.2, 2.0]]),
# ... 0.00001)
# >>> assert mask.sum() == 1 # one point in one triangle
>>> mesh = tfields.Mesh3D.grid((0, 1, 3), (1, 2, 3), (2, 3, 3))
>>> _ = mesh.tree
>>> assert hasattr(mesh, '_cache')
>>> assert 'mesh_tree' in mesh._cache
>>> mask = mesh.in_faces(tfields.Points3D([[0.2, 1.2, 2.0]]),
... 0.00001)
>>> assert mask.sum() == 1 # one point in one triangle
"""
raise ValueError("Broken feature. We are working on it!")
if not hasattr(self, "_cache"):
self._cache = {}
key = "mesh_tree"
......@@ -638,7 +633,7 @@ class Mesh3D(tfields.TensorMaps):
return_point_face_assignment=False,
):
"""
project the points of the tensor_field to a copy of the mesh
Project the points of the tensor_field to a copy of the mesh
and set the face values accord to the field to the maps field.
If no field is present in tensor_field, the number of points in a mesh
is counted.
......@@ -661,14 +656,15 @@ class Mesh3D(tfields.TensorMaps):
... [1, 2, 3, 4])
>>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,2,0], [1,2,0]],
... maps=[mp])
Projecting points onto the mesh gives the count
>>> points = tfields.Tensors([[0.5, 0.2, 0.0],
... [0.5, 0.02, 0.0],
... [0.5, 0.8, 0.0],
... [0.5, 0.8, 0.1]]) # not contained
Projecting points onto the mesh gives the count
>>> m_points = m.project(points, delta=0.01)
>>> assert m_points.maps[3].fields[0].equal([2, 1, 0, 0])
>>> list(m_points.maps[3].fields[0])
[2, 1, 0, 0]
TensorFields with arbitrary size are projected,
combinging the fields automatically
......@@ -680,8 +676,7 @@ class Mesh3D(tfields.TensorMaps):
... [[9, 9]] * 2])]
>>> tf = tfields.TensorFields(points, *fields)
>>> m_tf = m.project(tf, delta=0.01)
>>> assert m_tf.maps[3].fields[0].equal([2, 42, np.nan, np.nan],
... equal_nan=True)
>>> assert m_tf.maps[3].fields[0].equal([2, 42, np.nan, np.nan], equal_nan=True)
>>> assert m_tf.maps[3].fields[1].equal([[1, 2, 3],
... [3, 4, 5],
... [np.nan] * 3,
......@@ -730,18 +725,9 @@ class Mesh3D(tfields.TensorMaps):
if point_face_assignment is not None:
if len(point_face_assignment) != len(tensor_field):
raise ValueError("Template needs same lenght as tensor_field")
point_face_assignment = point_face_assignment
else:
mask = self.in_faces(tensor_field, delta=delta)
point_face_assignment = self.in_faces(tensor_field, delta=delta)
face_indices = np.arange(n_faces)
point_face_assignment = [
face_indices[mask[p_index, :]] for p_index in point_indices
]
point_face_assignment = np.array(
[fa if len(fa) != 0 else [-1] for fa in point_face_assignment]
)
point_face_assignment = point_face_assignment.reshape(-1)
point_face_assignment_set = set(point_face_assignment)
# merge the fields according to point_face_assignment
......
......@@ -7,6 +7,7 @@ Mail: daniel.boeckenhoff@ipp.mpg.de
part of tfields library
"""
import warnings
import typing
import numpy as np
import tfields
import logging
......@@ -529,54 +530,7 @@ class Triangles3D(tfields.TensorFields):
np.array: boolean mask, True where point in a triangle within delta
Examples:
>>> m = tfields.Mesh3D([[1,0,0], [0,1,0], [0,0,0]], faces=[[0, 1, 2]]);
>>> assert np.array_equal(
... m.triangles()._in_triangles(tfields.Points3D([[0.2, 0.2, 0]])),
... np.array([True], dtype=bool))
All Triangles are tested
>>> m2 = tfields.Mesh3D([[1,0,0], [0,1,0], [0,0,0], [4,0,0], [4, 4, 0], [8, 0, 0]],
... faces=[[0, 1, 2], [3, 4, 5]]);
>>> assert np.array_equal(
... m2.triangles()._in_triangles(np.array([0.2, 0.2, 0])),
... np.array([True, False], dtype=bool))
>>> assert np.array_equal(
... m2.triangles()._in_triangles(np.array([5, 2, 0])),
... np.array([False, True], dtype=bool))
delta allows to accept points that lie within delta orthogonal to the tringle plain
>>> assert np.array_equal(
... m2.triangles()._in_triangles(np.array([0.2, 0.2, 9000]), 0.0),
... np.array([False, False], dtype=bool))
>>> assert np.array_equal(
... m2.triangles()._in_triangles(np.array([0.2, 0.2, 0.1]), 0.2),
... np.array([ True, False], dtype=bool))
if you set delta to None, the minimal distance point(s) are accepted
>>> assert np.array_equal(
... m2.triangles()._in_triangles(np.array([0.2, 0.2, 0.1]), None),
... np.array([ True, False], dtype=bool))
If you define triangles that have colinear side vectors or in general lead to
not invertable matrices the you will always get False
>>> m3 = tfields.Mesh3D([[0,0,0], [2,0,0], [4,0,0], [0,1,0]],
... faces=[[0, 1, 2], [0, 1, 3]]);
>>> mask = m3.triangles()._in_triangles(np.array([0.2, 0.2, 0]), delta=0.3)
>>> assert np.array_equal(mask,
... np.array([False, True], dtype=bool))
Raises:
Wrong format of point will raise a ValueError
>>> m.triangles()._in_triangles(tfields.Points3D([[0.2, 0.2, 0], [0.2, 0.2, 0]]))
Traceback (most recent call last):
...
ValueError: point must be castable to shape 3 but is of shape (2, 3)
see tests/test_triangles_3d.py
"""
if self.ntriangles() == 0:
......@@ -628,7 +582,12 @@ class Triangles3D(tfields.TensorFields):
return np.array(barycentric_bools & orthogonal_acceptance)
def in_triangles(self, tensors, delta=0.0, assign_multiple=False):
def in_triangles(
self,
tensors,
delta: typing.Optional[float] = 0.0,
assign_multiple: bool = False,
) -> typing.Union[typing.List[typing.List[int]], np.array]:
"""
Barycentric method to obtain, which tensors are containes in any of the triangles
......@@ -636,61 +595,38 @@ class Triangles3D(tfields.TensorFields):
tensors (Points3D instance)
optional:
delta (:obj:`float` | :obj:`None`, optional):
:obj:`float`: normal distance to a triangle, that the points
delta:
:obj:`float`: Normal distance to a triangle, that the points
are concidered to be contained in the triangle.
:obj:`None`: find the minimum distance
default is 0.
assign_multiple (bool): if True, one point may belong to multiple
:obj:`None`: Find the minimum distance.
Default is 0.
assign_multiple: If True, one point may belong to multiple
triangles at the same time. In the other case the first
occurence will be True the other False
Returns:
np.ndarray: 2-d mask which is True, where a point is in a triangle
triangle index (axis = 1)->
points index
(axis = 0) 1 0 0
| 1 0 0
V 0 0 1
if assign_multiple == False:
there is just one True for each points index
triangle indices can have multiple true values
For Example, if you want to know the number of points in one
face, just sum over all points:
>> tris.in_triangles(poits).sum(axis=0)[face_index]
list: [index(or indices if assign_multiple) of triangle for point in tensors]
"""
indices = np.full(tensors.shape[0], -1, dtype=int)
if self.ntriangles() == 0:
return np.empty((tensors.shape[0], 0), dtype=bool)
return indices
masks = np.zeros((len(tensors), self.ntriangles()), dtype=bool)
with tensors.tmp_transform(self.coord_sys):
for i, point in enumerate(iter(tensors)):
masks[i] = self._in_triangles(point, delta)
if len(masks) == 0:
# empty sequence
return masks
if not assign_multiple:
"""
index of first faces containing the point for each point. This gets the
first membership index always. ignoring multiple triangle memberships
"""
faceMembershipIndices = np.argmax(masks, axis=1)
# True if point lies in any triangle
membershipBools = masks.sum(axis=1) != 0
masks = np.full(masks.shape, False, dtype=bool)
for j, valid in enumerate(membershipBools):
if valid:
masks[j, faceMembershipIndices[j]] = True
"""
masks is now the matrix as described in __doc__
"""
return masks
mask = self._in_triangles(point, delta)
if np.any(mask):
if assign_multiple:
index = np.argwhere(mask == np.amax(mask))
index.flatten().tolist()
indices.append(index)
else:
indices[i] = np.argmax(mask)
else:
if assign_multiple:
indices.append([-1])
else:
indices[i] = -1
return indices
def _on_edges(self, point):
"""
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment