diff --git a/tfields/__init__.py b/tfields/__init__.py
index 05f0dc305d875bc68a1ea7bc736eefad8302e546..3917788a337df04f437bb53842301529b3779f0e 100644
--- a/tfields/__init__.py
+++ b/tfields/__init__.py
@@ -1,3 +1,4 @@
+# flake8: noqa: F401
"""Top-level package of tfields."""
__author__ = """Daniel Böckenhoff"""
@@ -5,13 +6,13 @@ __email__ = "dboe@ipp.mpg.de"
__version__ = "0.3.3"
# methods:
-from tfields.core import dim, rank # NOQA
-from tfields.mask import evalf # NOQA
+from tfields.core import dim, rank
+from tfields.mask import evalf
# classes:
-from tfields.core import Tensors, TensorFields, TensorMaps, Container, Maps # NOQA
-from tfields.points3D import Points3D # NOQA
-from tfields.mesh3D import Mesh3D # NOQA
-from tfields.triangles3D import Triangles3D # NOQA
-from tfields.planes3D import Planes3D # NOQA
-from tfields.bounding_box import Node # NOQA
+from tfields.core import Tensors, TensorFields, TensorMaps, Container, Maps
+from tfields.points_3d import Points3D
+from tfields.mesh_3d import Mesh3D
+from tfields.triangles_3d import Triangles3D
+from tfields.planes_3d import Planes3D
+from tfields.bounding_box import Node
diff --git a/tfields/mesh3D.py b/tfields/mesh_3d.py
similarity index 98%
rename from tfields/mesh3D.py
rename to tfields/mesh_3d.py
index 2ed7f51ebc80e1253cc517179c03a54d5f1b1ec4..5965599018bd92496e1d583bf4a90428af6a9e98 100644
--- a/tfields/mesh3D.py
+++ b/tfields/mesh_3d.py
@@ -4,48 +4,48 @@
Author: Daniel Boeckenhoff
Mail: daniel.boeckenhoff@ipp.mpg.de
-part of tfields library
+Triangulated mesh class and methods
"""
+import logging
+import os
import numpy as np
import rna
import tfields
# obj imports
from tfields.lib.decorators import cached_property
-import logging
-import os
def _dist_from_plane(point, plane):
return plane["normal"].dot(point) + plane["d"]
-def _segment_plane_intersection(p0, p1, plane):
+def _segment_plane_intersection(point0, point1, plane):
"""
- Get the intersection between the ray spanned by p0 and p1 with the plane.
+ Get the intersection between the ray spanned by point0 and point1 with the plane.
Args:
- p0: array of length 3
- p1: array of length 3
+ point0: array of length 3
+ point1: array of length 3
plane:
Returns:
points, direction
points (list of arrays of length 3): 3d points
"""
- distance0 = _dist_from_plane(p0, plane)
- distance1 = _dist_from_plane(p1, plane)
- p0_on_plane = abs(distance0) < np.finfo(float).eps
- p1_on_plane = abs(distance1) < np.finfo(float).eps
+ distance0 = _dist_from_plane(point0, plane)
+ distance1 = _dist_from_plane(point1, plane)
+ point0_on_plane = abs(distance0) < np.finfo(float).eps
+ point1_on_plane = abs(distance1) < np.finfo(float).eps
points = []
direction = 0
- if p0_on_plane:
- points.append(p0)
+ if point0_on_plane:
+ points.append(point0)
- if p1_on_plane:
- points.append(p1)
+ if point1_on_plane:
+ points.append(point1)
# remove duplicate points
if len(points) > 1:
points = np.unique(points, axis=0)
- if p0_on_plane and p1_on_plane:
+ if point0_on_plane and point1_on_plane:
return points, direction
if distance0 * distance1 > np.finfo(float).eps:
@@ -54,14 +54,14 @@ def _segment_plane_intersection(p0, p1, plane):
direction = np.sign(distance0)
if abs(distance0) < np.finfo(float).eps:
return points, direction
- elif abs(distance1) < np.finfo(float).eps:
+ if abs(distance1) < np.finfo(float).eps:
return points, direction
if abs(distance0 - distance1) > np.finfo(float).eps:
t = distance0 / (distance0 - distance1)
else:
return points, direction
- points.append(p0 + t * (p1 - p0))
+ points.append(point0 + t * (point1 - point0))
# remove duplicate points
if len(points) > 1:
points = np.unique(points, axis=0)
diff --git a/tfields/planes3D.py b/tfields/planes_3d.py
similarity index 74%
rename from tfields/planes3D.py
rename to tfields/planes_3d.py
index ab94474ca8dbfab17d319e588bfa787a2194c7fa..00b705062df0a1e4730dd7ee0c71de77b8edfff0 100644
--- a/tfields/planes3D.py
+++ b/tfields/planes_3d.py
@@ -31,8 +31,10 @@ class Planes3D(tfields.TensorFields):
Returns:
list: list with sympy.Plane objects
"""
- return [sympy.Plane(point, normal_vector=vector)
- for point, vector in zip(self, self.fields[0])]
+ return [
+ sympy.Plane(point, normal_vector=vector)
+ for point, vector in zip(self, self.fields[0])
+ ]
def plot(self, **kwargs): # pragma: no cover
"""
@@ -42,14 +44,11 @@ class Planes3D(tfields.TensorFields):
centers = np.array(self)
norms = np.array(self.fields[0])
for i in range(len(self)):
- artists.append(
- rna.plotting.plot_plane(
- centers[i],
- norms[i],
- **kwargs))
+ artists.append(rna.plotting.plot_plane(centers[i], norms[i], **kwargs))
return artists
-if __name__ == '__main__': # pragma: no cover
+if __name__ == "__main__": # pragma: no cover
import doctest
+
doctest.testmod()
diff --git a/tfields/points3D.py b/tfields/points_3d.py
similarity index 93%
rename from tfields/points3D.py
rename to tfields/points_3d.py
index 2ce321180b541626db5e5414913ee7bea3868f58..e4f1ed412f2d2fc62ab1c98031870c9b16c2dc7c 100644
--- a/tfields/points3D.py
+++ b/tfields/points_3d.py
@@ -135,9 +135,10 @@ class Points3D(tfields.Tensors):
>>> assert p.equal(p_cart, atol=1e-15)
"""
+
def __new__(cls, tensors, **kwargs):
if not issubclass(type(tensors), Points3D):
- kwargs['dim'] = 3
+ kwargs["dim"] = 3
return super(Points3D, cls).__new__(cls, tensors, **kwargs)
def balls(self, radius, spacing=(5, 3)):
@@ -149,13 +150,15 @@ class Points3D(tfields.Tensors):
tfields.Mesh3D: Builds a sphere around each point with a resolution
defined by spacing and given radius
"""
- sphere = tfields.Mesh3D.grid((radius, radius, 1),
- (-np.pi, np.pi, spacing[0]),
- (-np.pi / 2, np.pi / 2, spacing[1]),
- coord_sys='spherical')
- sphere.transform('cartesian')
+ sphere = tfields.Mesh3D.grid(
+ (radius, radius, 1),
+ (-np.pi, np.pi, spacing[0]),
+ (-np.pi / 2, np.pi / 2, spacing[1]),
+ coord_sys="spherical",
+ )
+ sphere.transform("cartesian")
balls = []
- with self.tmp_transform('cartesian'):
+ with self.tmp_transform("cartesian"):
for point in self:
ball = sphere.copy()
ball += point
@@ -163,7 +166,8 @@ class Points3D(tfields.Tensors):
return tfields.Mesh3D.merged(*balls)
-if __name__ == '__main__': # pragma: no cover
+if __name__ == "__main__": # pragma: no cover
import doctest
+
doctest.testmod()
# doctest.run_docstring_examples(Points3D, globals())
diff --git a/tfields/triangles3D.py b/tfields/triangles_3d.py
similarity index 84%
rename from tfields/triangles3D.py
rename to tfields/triangles_3d.py
index 3d0d15d33b34217cd3680082726277badfdb61c8..3483912f920092c8063ff3e19660d10e942e16a4 100644
--- a/tfields/triangles3D.py
+++ b/tfields/triangles_3d.py
@@ -39,14 +39,15 @@ class Triangles3D(tfields.TensorFields):
"""
def __new__(cls, tensors, *fields, **kwargs):
- kwargs['dim'] = 3
- kwargs['rigid'] = False
+ kwargs["dim"] = 3
+ kwargs["rigid"] = False
obj = super(Triangles3D, cls).__new__(cls, tensors, *fields, **kwargs)
if not len(obj) % 3 == 0:
- warnings.warn("Input object of size({0}) has no divider 3 and"
- " does not describe triangles."
- .format(len(obj)))
+ warnings.warn(
+ "Input object of size({0}) has no divider 3 and"
+ " does not describe triangles.".format(len(obj))
+ )
return obj
def ntriangles(self):
@@ -105,13 +106,16 @@ class Triangles3D(tfields.TensorFields):
# apply triangle index to fields
if tri_index is not None:
- item.fields = [field.__getitem__(tri_index)
- for field in item.fields]
+ item.fields = [
+ field.__getitem__(tri_index) for field in item.fields
+ ]
else:
item = tfields.Tensors(item)
except IndexError:
- logging.warning("Index error occured for field.__getitem__. Error "
- "message: {err}".format(**locals()))
+ logging.warning(
+ "Index error occured for field.__getitem__. Error "
+ "message: {err}".format(**locals())
+ )
return item
@classmethod
@@ -121,6 +125,7 @@ class Triangles3D(tfields.TensorFields):
Given a path to a stl file, construct the object
"""
import stl.mesh
+
triangles = stl.mesh.Mesh.from_file(path)
obj = cls(triangles.vectors.reshape(-1, 3))
return obj
@@ -128,12 +133,13 @@ class Triangles3D(tfields.TensorFields):
@classmethod
def merged(cls, *objects, **kwargs):
with warnings.catch_warnings():
- warnings.filterwarnings('ignore')
+ warnings.filterwarnings("ignore")
obj = super(Triangles3D, cls).merged(*objects, **kwargs)
if not len(obj) % 3 == 0:
- warnings.warn("Input object of size({0}) has no divider 3 and"
- " does not describe triangles."
- .format(len(obj)))
+ warnings.warn(
+ "Input object of size({0}) has no divider 3 and"
+ " does not describe triangles.".format(len(obj))
+ )
return obj
def evalf(self, expression=None, coord_sys=None):
@@ -189,9 +195,7 @@ class Triangles3D(tfields.TensorFields):
Returns:
tfields.Mesh3D
"""
- mp = tfields.TensorFields(
- np.arange(len(self)).reshape((-1, 3)),
- *self.fields)
+ mp = tfields.TensorFields(np.arange(len(self)).reshape((-1, 3)), *self.fields)
mesh = tfields.Mesh3D(self, maps=[mp])
return mesh.cleaned(stale=False) # stale vertices can not occure here
@@ -237,30 +241,33 @@ class Triangles3D(tfields.TensorFields):
# define 3 vectors building the triangle, transform it back and take their norm
if not np.array_equal(transform, np.eye(3)):
- a = np.linalg.norm(np.linalg.solve(transform.T,
- (self[aIndices, :] -
- self[bIndices, :]).T),
- axis=0)
- b = np.linalg.norm(np.linalg.solve(transform.T,
- (self[aIndices, :] -
- self[cIndices, :]).T),
- axis=0)
- c = np.linalg.norm(np.linalg.solve(transform.T,
- (self[bIndices, :] -
- self[cIndices, :]).T),
- axis=0)
+ a = np.linalg.norm(
+ np.linalg.solve(
+ transform.T, (self[aIndices, :] - self[bIndices, :]).T
+ ),
+ axis=0,
+ )
+ b = np.linalg.norm(
+ np.linalg.solve(
+ transform.T, (self[aIndices, :] - self[cIndices, :]).T
+ ),
+ axis=0,
+ )
+ c = np.linalg.norm(
+ np.linalg.solve(
+ transform.T, (self[bIndices, :] - self[cIndices, :]).T
+ ),
+ axis=0,
+ )
else:
- a = np.linalg.norm(self[aIndices, :] - self[bIndices, :],
- axis=1)
- b = np.linalg.norm(self[aIndices, :] - self[cIndices, :],
- axis=1)
- c = np.linalg.norm(self[bIndices, :] - self[cIndices, :],
- axis=1)
+ a = np.linalg.norm(self[aIndices, :] - self[bIndices, :], axis=1)
+ b = np.linalg.norm(self[aIndices, :] - self[cIndices, :], axis=1)
+ c = np.linalg.norm(self[bIndices, :] - self[cIndices, :], axis=1)
# sort by length for numerical stability
lengths = np.concatenate(
- (a.reshape(-1, 1), b.reshape(-1, 1), c.reshape(-1, 1)),
- axis=1)
+ (a.reshape(-1, 1), b.reshape(-1, 1), c.reshape(-1, 1)), axis=1
+ )
lengths.sort()
a, b, c = lengths.T
@@ -300,18 +307,22 @@ class Triangles3D(tfields.TensorFields):
"""
pointA, pointB, pointC = self.corners()
- a = np.linalg.norm(pointC - pointB, axis=1) # side length of side opposite to pointA
+ a = np.linalg.norm(
+ pointC - pointB, axis=1
+ ) # side length of side opposite to pointA
b = np.linalg.norm(pointC - pointA, axis=1)
c = np.linalg.norm(pointB - pointA, axis=1)
bary1 = a ** 2 * (b ** 2 + c ** 2 - a ** 2)
bary2 = b ** 2 * (a ** 2 + c ** 2 - b ** 2)
bary3 = c ** 2 * (a ** 2 + b ** 2 - c ** 2)
- matrices = np.concatenate((pointA, pointB, pointC), axis=1).reshape(pointA.shape + (3,))
+ matrices = np.concatenate((pointA, pointB, pointC), axis=1).reshape(
+ pointA.shape + (3,)
+ )
# transpose the inner matrix
- matrices = np.einsum('...ji', matrices)
+ matrices = np.einsum("...ji", matrices)
vectors = np.array((bary1, bary2, bary3)).T
# matrix vector product for matrices and vectors
- P = np.einsum('...ji,...i', matrices, vectors)
+ P = np.einsum("...ji,...i", matrices, vectors)
P /= vectors.sum(axis=1).reshape((len(vectors), 1))
return tfields.Points3D(P)
@@ -324,7 +335,9 @@ class Triangles3D(tfields.TensorFields):
nT = self.ntriangles()
mat = np.ones((1, 3)) / 3.0
# matrix product calculatesq center of all triangles
- return tfields.Points3D(np.dot(mat, self.reshape(nT, 3, 3))[0], coord_sys=self.coord_sys)
+ return tfields.Points3D(
+ np.dot(mat, self.reshape(nT, 3, 3))[0], coord_sys=self.coord_sys
+ )
"""
Old version:
@@ -387,7 +400,7 @@ class Triangles3D(tfields.TensorFields):
np.place(norms, norms == 0.0, 1.0)
return vectors / norms
- def _baricentric(self, point, delta=0.):
+ def _baricentric(self, point, delta=0.0):
"""
Determine baricentric coordinates like
@@ -446,14 +459,22 @@ class Triangles3D(tfields.TensorFields):
ap = point - a
# matrix vector product for matrices and vectors
- barCoords = np.einsum('...ji,...i', self._baricentric_matrix, ap)
+ barCoords = np.einsum("...ji,...i", self._baricentric_matrix, ap)
with warnings.catch_warnings():
# python2.7
- warnings.filterwarnings('ignore', message="invalid value encountered in divide")
- warnings.filterwarnings('ignore', message="divide by zero encountered in divide")
+ warnings.filterwarnings(
+ "ignore", message="invalid value encountered in divide"
+ )
+ warnings.filterwarnings(
+ "ignore", message="divide by zero encountered in divide"
+ )
# python3.x
- warnings.filterwarnings('ignore', message="invalid value encountered in true_divide")
- warnings.filterwarnings('ignore', message="divide by zero encountered in true_divide")
+ warnings.filterwarnings(
+ "ignore", message="invalid value encountered in true_divide"
+ )
+ warnings.filterwarnings(
+ "ignore", message="divide by zero encountered in true_divide"
+ )
barCoords[:, 2] /= delta # allow devide by 0.
return barCoords
@@ -468,15 +489,16 @@ class Triangles3D(tfields.TensorFields):
norm = np.cross(ab, ac)
normLen = np.linalg.norm(norm, axis=1)
normLen = normLen.reshape((1,) + normLen.shape)
- colinear_mask = (normLen == 0)
+ colinear_mask = normLen == 0
with warnings.catch_warnings():
- warnings.filterwarnings('ignore', category=np.VisibleDeprecationWarning)
+ warnings.filterwarnings("ignore", category=np.VisibleDeprecationWarning)
# prevent divide by 0
- norm[np.where(~colinear_mask.T)] = \
+ norm[np.where(~colinear_mask.T)] = (
norm[np.where(~colinear_mask.T)] / normLen.T[np.where(~colinear_mask.T)]
+ )
matrix = np.concatenate((ab, ac, norm), axis=1).reshape(ab.shape + (3,))
- matrix = np.einsum('...ji', matrix) # transpose the inner matrix
+ matrix = np.einsum("...ji", matrix) # transpose the inner matrix
# invert matrix if possible
# matrixI = np.linalg.inv(matrix) # one line variant without exception
@@ -486,13 +508,14 @@ class Triangles3D(tfields.TensorFields):
matrixI.append(np.linalg.inv(mat))
except np.linalg.linalg.LinAlgError as e:
if str(e) == "Singular matrix":
- warnings.warn("Singular matrix: Could not invert matrix "
- "{0}. Return nan matrix."
- .format(str(mat).replace('\n', ',')))
+ warnings.warn(
+ "Singular matrix: Could not invert matrix "
+ "{0}. Return nan matrix.".format(str(mat).replace("\n", ","))
+ )
matrixI.append(np.full((3, 3), np.nan))
return np.array(matrixI)
- def _in_triangles(self, point, delta=0.):
+ def _in_triangles(self, point, delta=0.0):
"""
Barycentric method to optain, wheter a point is in any of the triangles
@@ -562,46 +585,50 @@ class Triangles3D(tfields.TensorFields):
try:
point = np.reshape(point, 3)
except ValueError:
- raise ValueError("point must be castable to shape 3 but is of shape {0}"
- .format(point.shape))
+ raise ValueError(
+ "point must be castable to shape 3 but is of shape {0}".format(
+ point.shape
+ )
+ )
# min_dist_method switch if delta is None
if delta is None:
- delta = 1.
+ delta = 1.0
min_dist_method = True
else:
min_dist_method = False
u, v, w = self._baricentric(point, delta=delta).T
- if delta == 0.:
- w[np.isnan(w)] = 0. # division by 0 in baricentric makes w = 0 nan.
+ if delta == 0.0:
+ w[np.isnan(w)] = 0.0 # division by 0 in baricentric makes w = 0 nan.
with warnings.catch_warnings():
warnings.filterwarnings(
- 'ignore',
- message="invalid value encountered in less_equal")
- barycentric_bools = ((v <= 1.) & (v >= 0.))\
- & ((u <= 1.) & (u >= 0.))\
- & ((v + u <= 1.0))
+ "ignore", message="invalid value encountered in less_equal"
+ )
+ barycentric_bools = (
+ ((v <= 1.0) & (v >= 0.0)) & ((u <= 1.0) & (u >= 0.0)) & ((v + u <= 1.0))
+ )
if all(~barycentric_bools):
return barycentric_bools
if min_dist_method:
orthogonal_acceptance = np.full(
- barycentric_bools.shape, False, dtype=bool)
+ barycentric_bools.shape, False, dtype=bool
+ )
closest_indices = np.argmin(abs(w)[barycentric_bools])
# transform the indices to the whole array, not only the
# barycentric_bools selection
- closest_indices = np.arange(
- len(barycentric_bools)
- )[barycentric_bools][closest_indices]
+ closest_indices = np.arange(len(barycentric_bools))[barycentric_bools][
+ closest_indices
+ ]
orthogonal_acceptance[closest_indices] = True
else:
- orthogonal_acceptance = (abs(w) <= 1)
+ orthogonal_acceptance = abs(w) <= 1
return np.array(barycentric_bools & orthogonal_acceptance)
- def in_triangles(self, tensors, delta=0., assign_multiple=False):
+ def in_triangles(self, tensors, delta=0.0, assign_multiple=False):
"""
Barycentric method to obtain, which tensors are containes in any of the triangles
@@ -654,7 +681,7 @@ class Triangles3D(tfields.TensorFields):
"""
faceMembershipIndices = np.argmax(masks, axis=1)
# True if point lies in any triangle
- membershipBools = (masks.sum(axis=1) != 0)
+ membershipBools = masks.sum(axis=1) != 0
masks = np.full(masks.shape, False, dtype=bool)
for j, valid in enumerate(membershipBools):
@@ -693,12 +720,14 @@ class Triangles3D(tfields.TensorFields):
... np.array([False, False, True], dtype=bool))
"""
- u, v, w = self._baricentric(point, 1.).T
+ u, v, w = self._baricentric(point, 1.0).T
- orthogonal_acceptance = (w == 0) # point should lie in triangle
- barycentric_bools = (((0. <= v) & (v <= 1.)) & (u == 0.)) | \
- (((0. <= u) & (u <= 1.)) & (v == 0.)) | \
- (v + u == 1.0)
+ orthogonal_acceptance = w == 0 # point should lie in triangle
+ barycentric_bools = (
+ (((0.0 <= v) & (v <= 1.0)) & (u == 0.0))
+ | (((0.0 <= u) & (u <= 1.0)) & (v == 0.0))
+ | (v + u == 1.0)
+ )
return np.array(barycentric_bools & orthogonal_acceptance)
def _weights(self, weights, rigid=False):
@@ -722,6 +751,7 @@ class Triangles3D(tfields.TensorFields):
return super(Triangles3D, self)._weights(weights, rigid=rigid)
-if __name__ == '__main__': # pragma: no cover
+if __name__ == "__main__": # pragma: no cover
import doctest
+
doctest.testmod()