Commit a6984734 authored by dboe's avatar dboe
Browse files

renaming 3D -> _3d

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