From 9469dd0ecb66096fd98bdfdf348dac37e1e5105c Mon Sep 17 00:00:00 2001
From: dboe <dboe@ipp.mpg.de>
Date: Wed, 7 Jul 2021 17:24:21 +0200
Subject: [PATCH] added test for _in_triangles and in_triangles
---
tests/test_triangles_3d.py | 108 ++++++++++++++++++++++++++++++++
tfields/bounding_box.py | 18 +++---
tfields/mesh_3d.py | 54 ++++++----------
tfields/triangles_3d.py | 124 +++++++++----------------------------
4 files changed, 165 insertions(+), 139 deletions(-)
create mode 100644 tests/test_triangles_3d.py
diff --git a/tests/test_triangles_3d.py b/tests/test_triangles_3d.py
new file mode 100644
index 0000000..e0e2c4f
--- /dev/null
+++ b/tests/test_triangles_3d.py
@@ -0,0 +1,108 @@
+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
diff --git a/tfields/bounding_box.py b/tfields/bounding_box.py
index b3d0145..989c841 100644
--- a/tfields/bounding_box.py
+++ b/tfields/bounding_box.py
@@ -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__":
diff --git a/tfields/mesh_3d.py b/tfields/mesh_3d.py
index afcebff..0c8d68e 100644
--- a/tfields/mesh_3d.py
+++ b/tfields/mesh_3d.py
@@ -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
diff --git a/tfields/triangles_3d.py b/tfields/triangles_3d.py
index 3483912..2d263e5 100644
--- a/tfields/triangles_3d.py
+++ b/tfields/triangles_3d.py
@@ -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):
"""
--
GitLab