From 64f4270c4c5556b012c051736131be05f317e41f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20B=C3=B6ckenhoff=20=28Laptop=29?= <dboe@ipp.mpg.de> Date: Wed, 13 Jun 2018 18:17:32 +0200 Subject: [PATCH] mesh3D runs; triangles mostly --- tfields/__init__.py | 4 +- tfields/core.py | 124 ++++++-- tfields/mask.py | 15 +- tfields/mesh3D.py | 52 ++-- tfields/triangles3D.py | 670 +++++++++++++---------------------------- 5 files changed, 339 insertions(+), 526 deletions(-) diff --git a/tfields/__init__.py b/tfields/__init__.py index 473582f..4a5ef44 100644 --- a/tfields/__init__.py +++ b/tfields/__init__.py @@ -6,10 +6,10 @@ from .lib import * # __all__ = ['core', 'points3D'] from .core import Tensors, TensorFields, TensorMaps from .points3D import Points3D -from .mask import getMask +from .mask import evalf # methods: -from .mask import getMask # NOQA +from .mask import evalf # NOQA from .lib import * # NOQA # classes: diff --git a/tfields/core.py b/tfields/core.py index 156d0ff..e0a6141 100644 --- a/tfields/core.py +++ b/tfields/core.py @@ -692,7 +692,7 @@ class Tensors(AbstractNdarray): if condition is None: condition = np.array([True for i in range(len(self))]) elif isinstance(condition, sympy.Basic): - condition = self.getMask(condition) + condition = self.evalf(condition) if isinstance(coordinate, list) or isinstance(coordinate, tuple): for c in coordinate: self.mirror(c, condition) @@ -803,7 +803,7 @@ class Tensors(AbstractNdarray): raise ValueError("Multiple occurences of value {}" .format(tensor)) - def getMoment(self, moment): + def moments(self, moment): """ Returns: Moments of the distribution. @@ -813,9 +813,9 @@ class Tensors(AbstractNdarray): Args: moment (int): n-th moment """ - return tfields.lib.stats.getMoment(self, moment) + return tfields.lib.stats.moments(self, moment) - def closestPoints(self, other, **kwargs): + def closest(self, other, **kwargs): """ Args: other (Tensors): closest points to what? -> other @@ -827,26 +827,27 @@ class Tensors(AbstractNdarray): >>> m = tfields.Tensors([[1,0,0], [0,1,0], [1,1,0], [0,0,1], ... [1,0,1]]) >>> p = tfields.Tensors([[1.1,1,0], [0,0.1,1], [1,0,1.1]]) - >>> p.closestPoints(m) + >>> p.closest(m) array([2, 3, 4]) """ with other.tmp_transform(self.coordSys): # balanced_tree option gives huge speedup! - kdTree = sp.spatial.cKDTree(other, 1000, - balanced_tree=False) - res = kdTree.query(self, **kwargs) + kd_tree = sp.spatial.cKDTree(other, 1000, + balanced_tree=False) + res = kd_tree.query(self, **kwargs) array = res[1] return array - def getMask(self, cutExpression=None, coordSys=None): + def evalf(self, expression=None, coordSys=None): """ Args: - cutExpression (sympy logical expression) - coordSys (str): coordSys to evaluate the expression in. - Returns: np.array of dtype bool with lenght of number of points in self. - This array is True, where cutExpression evaluates True. + expression (sympy logical expression) + coordSys (str): coordSys to evalfuate the expression in. + Returns: + np.ndarray: mask of dtype bool with lenght of number of points in self. + This array is True, where expression evalfuates True. Examples: >>> import tfields >>> import numpy @@ -854,33 +855,33 @@ class Tensors(AbstractNdarray): >>> x, y, z = sympy.symbols('x y z') >>> p = tfields.Points3D([[1., 2., 3.], [4., 5., 6.], [1, 2, -6], ... [-5, -5, -5], [1,0,-1], [0,1,-1]]) - >>> np.array_equal(p.getMask(x > 0), + >>> np.array_equal(p.evalf(x > 0), ... [True, True, True, False, True, False]) True And combination - >>> np.array_equal(p.getMask((x > 0) & (y < 3)), + >>> np.array_equal(p.evalf((x > 0) & (y < 3)), ... [True, False, True, False, True, False]) True Or combination - >>> np.array_equal(p.getMask((x > 0) | (y > 3)), + >>> np.array_equal(p.evalf((x > 0) | (y > 3)), ... [True, True, True, False, True, False]) True """ coords = sympy.symbols('x y z') with self.tmp_transform(coordSys or self.coordSys): - mask = tfields.getMask(self, cutExpression, coords=coords) + mask = tfields.evalf(self, expression, coords=coords) return mask - def cut(self, cutExpression, coordSys=None): + def cut(self, expression, coordSys=None): """ Default cut method for Points3D. Works on a copy. Args: - cutExpression (sympy logical expression): logical expression which will be evaluated. + expression (sympy logical expression): logical expression which will be evalfuated. use symbols x, y and z - coordSys (str): coordSys to evaluate the expression in. + coordSys (str): coordSys to evalfuate the expression in. Examples: >>> import tfields >>> import sympy @@ -903,7 +904,7 @@ class Tensors(AbstractNdarray): """ if len(self) == 0: return self.copy() - mask = self.getMask(cutExpression, coordSys=coordSys or self.coordSys) + mask = self.evalf(expression, coordSys=coordSys or self.coordSys) mask.astype(bool) inst = self[mask].copy() return inst @@ -933,15 +934,18 @@ class Tensors(AbstractNdarray): other.transform(self.coordSys) return sp.spatial.distance.cdist(self, other, **kwargs) - def minDistances(self, other=None, **kwargs): + def min_dists(self, other=None, **kwargs): """ Args: - other(array or None) + other(array | None): if None: closest distance to self **kwargs: memory_saving (bool): for very large array comparisons default False ... rest is forwarded to sp.spatial.distance.cdist + Returns: + np.array: minimal distances of self to other + Examples: >>> import tfields @@ -950,13 +954,13 @@ class Tensors(AbstractNdarray): ... (0, 2, 3), ... (0, 0, 1)) >>> p[4,2] = 1 - >>> dMin = p.minDistances() + >>> dMin = p.min_dists() >>> expected = [1] * 9 >>> expected[4] = np.sqrt(2) >>> np.array_equal(dMin, expected) True - >>> dMin2 = p.minDistances(memory_saving=True) + >>> dMin2 = p.min_dists(memory_saving=True) >>> bool((dMin2 == dMin).all()) True @@ -1001,6 +1005,60 @@ class Tensors(AbstractNdarray): distsInEpsilon = dists <= epsilon return [indices[die] for die in distsInEpsilon] + def _weights(self, weights, rigid=True): + """ + transformer method for weights inputs. + Args: + weights (np.ndarray | None): + If weights is None, use np.ones + Otherwise just pass the weights. + rigid (bool): demand equal weights and tensor length + Returns: + weight array + """ + # set weights to 1.0 if weights is None + if weights is None: + weights = np.ones(len(self)) + if rigid: + if not len(weights) == len(self): + raise ValueError("Equal number of weights as tensors demanded.") + return weights + + def cov_eig(self, weights=None): + """ + Calculate the covariance eigenvectors with lenghts of eigenvalues + Args: + weights (np.array | int | None): index to scalars to weight with + """ + # weights = self.getNormedWeightedAreas(weights=weights) + weights = self._weights(weights) + cov = np.cov(self.T, + ddof=0, + aweights=weights) + # calculate eigenvalues and eigenvectors of covariance + evalfs, evecs = np.linalg.eigh(cov) + idx = evalfs.argsort()[::-1] + evalfs = evalfs[idx] + evecs = evecs[:, idx] + e = np.concatenate((evecs, evalfs.reshape(1, 3))) + return e.T.reshape(12, ) + + def main_axes(self, weights=None): + """ + Returns: + Main Axes eigen-vectors + """ + # weights = self.getNormedWeightedAreas(weights=weights) + weights = self._weights(weights) + mean = self.moments(1) + relative_coords = self - mean + cov = np.cov(relative_coords.T, + ddof=0, + aweights=weights) + # calculate eigenvalues and eigenvectors of covariance + evalfs, evecs = np.linalg.eigh(cov) + return (evecs * evalfs.T).T + class TensorFields(Tensors): """ @@ -1068,7 +1126,7 @@ class TensorFields(Tensors): This error can be suppressed by setting rigid=False >>> loose = tfields.TensorFields([1, 2, 3], [3], rigid=False) - + >>> assert len(loose) != 1 """ __slots__ = ['coordSys', 'fields'] @@ -1192,6 +1250,18 @@ class TensorFields(Tensors): mask &= field.equal(other.fields[i], **kwargs) return mask + def _weights(self, weights, rigid=True): + """ + Expansion of Tensors._weights with integer inputs + Args: + weights (np.ndarray | int | None): + if weights is int: use field at index <weights> + else: see Tensors._weights + """ + if isinstance(weights, int) : + weights = self.fields[weights] + return super(TensorFields, self)._weights(weights, rigid=rigid) + class TensorMaps(TensorFields): """ @@ -1431,7 +1501,7 @@ class TensorMaps(TensorFields): >>> m = TensorMaps([[1,2,3], [3,3,3], [0,0,0], [5,6,7]], ... maps=[[[0, 1, 2], [1, 2, 3]]]) >>> from sympy.abc import x,y,z - >>> vertexMask = m.getMask(z < 6) + >>> vertexMask = m.evalf(z < 6) >>> faceMask = m.to_maps_masks(vertexMask) >>> assert np.array_equal(faceMask, [[True, False]]) diff --git a/tfields/mask.py b/tfields/mask.py index 14b85f7..62a7389 100644 --- a/tfields/mask.py +++ b/tfields/mask.py @@ -11,15 +11,16 @@ import numpy as np import sympy -def getMask(array, cutExpression=None, coords=None): +def evalf(array, cutExpression=None, coords=None): """ Linking sympy and numpy by retrieving a mask according to the cutExpression Args: array (numpy ndarray) cutExpression (sympy logical expression) - coordSys (str): coordSys to evaluate the expression in. - Returns: np.array which is True, where cutExpression evaluates True. + coordSys (str): coordSys to evalfuate the expression in. + Returns: + np.array: mask which is True, where cutExpression evalfuates True. Examples: >>> import sympy >>> import numpy as np @@ -27,20 +28,20 @@ def getMask(array, cutExpression=None, coords=None): >>> x, y, z = sympy.symbols('x y z') >>> a = np.array([[1., 2., 3.], [4., 5., 6.], [1, 2, -6], [-5, -5, -5], [1,0,-1], [0,1,-1]]) - >>> assert np.array_equal(tfields.getMask(a, x > 0), + >>> assert np.array_equal(tfields.evalf(a, x > 0), ... np.array([ True, True, True, False, True, False])) And combination - >>> assert np.array_equal(tfields.getMask(a, (x > 0) & (y < 3)), + >>> assert np.array_equal(tfields.evalf(a, (x > 0) & (y < 3)), ... np.array([True, False, True, False, True, False])) Or combination - >>> assert np.array_equal(tfields.getMask(a, (x > 0) | (y > 3)), + >>> assert np.array_equal(tfields.evalf(a, (x > 0) | (y > 3)), ... np.array([True, True, True, False, True, False])) If array of other shape than (?, 3) is given, the coords need to be specified >>> a0, a1 = sympy.symbols('a0 a1') - >>> assert np.array_equal(tfields.getMask([[0., 1.], [-1, 3]], a1 > 2, + >>> assert np.array_equal(tfields.evalf([[0., 1.], [-1, 3]], a1 > 2, ... coords=[a0, a1]), ... np.array([False, True], dtype=bool)) diff --git a/tfields/mesh3D.py b/tfields/mesh3D.py index aea8dd1..6326811 100644 --- a/tfields/mesh3D.py +++ b/tfields/mesh3D.py @@ -36,7 +36,7 @@ def fields_to_scalars(fields): return np.array(fields) def faces_to_maps(faces, *fields): - return [tfields.TensorFields(faces, *fields, dtype=int)] + return [tfields.TensorFields(faces, *fields, dtype=int, dim=3)] def maps_to_faces(maps): if len(maps) == 0: @@ -74,10 +74,9 @@ class Mesh3D(tfields.TensorMaps): going from Mesh3D to Triangles3D instance is easy and will be cached. >>> m = tfields.Mesh3D([[1,0,0], [0,1,0], [0,0,0]], faces=[[0, 1, 2]]); - >>> m.triangles - Triangles3D([[ 1., 0., 0.], - [ 0., 1., 0.], - [ 0., 0., 0.]]) + >>> assert m.triangles.equal(tfields.Triangles3D([[ 1., 0., 0.], + ... [ 0., 1., 0.], + ... [ 0., 0., 0.]])) a list of scalars is assigned to each face >>> mScalar = tfields.Mesh3D([[1,0,0], [0,1,0], [0,0,0]], faces=[[0, 1, 2]], faceScalars=[.5]); @@ -95,9 +94,7 @@ class Mesh3D(tfields.TensorMaps): ... [ 2., 0., 0.], ... [ 0., 3., 0.]]) True - >>> msum.faces - array([[0, 1, 2], - [3, 4, 5]]) + >>> assert np.array_equal(msum.faces, [[0, 1, 2], [3, 4, 5]]) Saving and reading >>> from tempfile import NamedTemporaryFile @@ -112,9 +109,9 @@ class Mesh3D(tfields.TensorMaps): """ def __new__(cls, tensors, **kwargs): - fields = kwargs.pop('fields', []) if not issubclass(type(tensors), Mesh3D): kwargs['dim'] = 3 + fields = kwargs.pop('fields', []) faces = kwargs.pop('faces', None) faceScalars = kwargs.pop('faceScalars', []) maps = kwargs.pop('maps', None) @@ -288,7 +285,10 @@ class Mesh3D(tfields.TensorMaps): """ if self.faces.size == 0: return tfields.Triangles3D([]) - return tfields.Triangles3D(self[self.faces.flatten()], triangleScalars=self.faceScalars) + tris = tfields.Tensors.merged(*[self[mp.flatten()] for mp in self.maps]) + map_fields = [mp.fields for mp in self.maps] + fields = [tfields.Tensors.merged(*fields) for fields in zip(*map_fields)] + return tfields.Triangles3D(tris, *fields) @decoTools.cached_property() def planes(self): @@ -313,22 +313,22 @@ class Mesh3D(tfields.TensorMaps): def createFromTxtFile(cls, filePath): return tfields.Triangles3D.createFromTxtFile(filePath).getMesh3D() - def __getattr__(self, name): - """ - getter methods are forwarded to self.triangles - Examples: - >>> m = tfields.Mesh3D([]); - >>> m.getAreas - <bound method Triangles3D.getAreas of Triangles3D([], shape=(0, 3), dtype=float64)> - - """ - if name.startswith('get'): - if not hasattr(tfields.Triangles3D, name) or name == 'getMask': - raise AttributeError("Could not forward attribute {0}".format(name)) - else: - return getattr(self.triangles, name) - else: - raise AttributeError("No attribute with name {0}".format(name)) + # def __getattr__(self, name): + # """ + # getter methods are forwarded to self.triangles + # Examples: + # >>> m = tfields.Mesh3D([]); + # >>> m.getAreas + # <bound method Triangles3D.getAreas of Triangles3D([], shape=(0, 3), dtype=float64)> + + # """ + # if name.startswith('get'): + # if not hasattr(tfields.Triangles3D, name) or name == 'getMask': + # raise AttributeError("Could not forward attribute {0}".format(name)) + # else: + # return getattr(self.triangles, name) + # else: + # raise AttributeError("No attribute with name {0}".format(name)) def getNFaces(self): return self.faces.shape[0] diff --git a/tfields/triangles3D.py b/tfields/triangles3D.py index 1a97570..8ac2612 100644 --- a/tfields/triangles3D.py +++ b/tfields/triangles3D.py @@ -6,15 +6,11 @@ Mail: daniel.boeckenhoff@ipp.mpg.de part of tfields library """ -import numpy as np import warnings +import logging +import numpy as np import tfields -import ioTools import decoTools -import pyTools -import loggingTools - -logger = loggingTools.Logger(__name__) class Triangles3D(tfields.TensorFields): @@ -22,42 +18,37 @@ class Triangles3D(tfields.TensorFields): """ Points3D child restricted to n * 3 Points. 3 Points always group together to one triangle. Examples: + >>> import tfields + >>> import numpy as np >>> t = tfields.Triangles3D([[1,2,3], [3,3,3], [0,0,0]]) - >>> t.triangleScalars - array([], shape=(1, 0), dtype=float64) - >>> _ = t.transform(tfields.bases.CYLINDER) - >>> a = [tfields.Points3D([[1, 2, 3]], coordSys=tfields.bases.CYLINDER), - ... tfields.Points3D([[3] * 3]), - ... tfields.Points3D([[5, 1, 3]])] - >>> tfields.Triangles3D.merged(*a, coordSys=tfields.bases.CARTESIAN) - Triangles3D([[-0.41614684, 0.90929743, 3. ], - [ 3. , 3. , 3. ], - [ 5. , 1. , 3. ]]) - - You can add triangleScalars to each triangle - >>> t.appendScalars([1]) + You can add fields to each triangle + >>> t = tfields.Triangles3D(t, tfields.Tensors([42])) + >>> assert t.fields[0].equal([42]) """ - # __slots__ = ['coordSys', 'triangleScalars', '_cache'] - def __new__(cls, tensors, *fields, **kwargs): - if not issubclass(type(tensors), Triangles3D): - 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: - raise TypeError("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 @classmethod def merged(cls, *objects, **kwargs): - obj = tfields.TensorFields.merged(*objects, **kwargs) - return cls.__new__(cls, obj) + with warnings.catch_warnings(): + warnings.filterwarnings('ignore') + obj = super(Triangles3D, cls).merged(cls, *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))) + return obj @property def triangleScalars(self): @@ -70,48 +61,29 @@ class Triangles3D(tfields.TensorFields): def triangleScalars(self, scalars): self.fields = tfields.scalars_to_fields(scalars) - def saveTxt(self, filePath): - self.getMesh3D().save(filePath) - - def appendScalars(self, scalars): - """ - Examples: - >>> t = tfields.Triangles3D([[1,2,3], [3,3,3], [0,0,0]], [[1001]]) - >>> t.appendScalars([42]) - >>> t.triangleScalars - array([[1001, 42]]) - - """ - self.triangleScalars = np.concatenate([self.triangleScalars, np.array([scalars]).T], 1) - - def getScalarArrays(self): - """ - Returns: - np.array: transposed scalar arrays - """ - return self.triangleScalars.T - - def getNTriangles(self): + def get_ntriangles(self): """ Returns: int: number of triangles """ - return self.getNPoints() // 3 + return len(self) // 3 - def getMask(self, cutExpression, coordSys=None): + def evalf(self, expression=None, coordSys=None): """ Triangle3D implementation """ - mask = super(Triangles3D, self).getMask(cutExpression, coordSys=coordSys) - mask = mask.reshape((self.getNTriangles(), 3)) + mask = super(Triangles3D, self).evalf(expression, coordSys=coordSys) + mask = mask.reshape((self.get_ntriangles(), 3)) mask = mask.all(axis=1) - mask = np.array([mask] * 3).T.reshape((self.getNPoints())) + mask = np.array([mask] * 3).T.reshape((len(self))) return mask - def cut(self, cutExpression, coordSys=None): + def cut(self, expression, coordSys=None): """ Default cut method for Triangles3D Examples: + >>> import sympy + >>> import tfields >>> x, y, z = sympy.symbols('x y z') >>> t = tfields.Triangles3D([[1., 2., 3.], [-4., 5., 6.], [1, 2, -6], ... [5, -5, -5], [1,0,-1], [0,1,-1], [-5, -5, -5], @@ -129,83 +101,81 @@ class Triangles3D(tfields.TensorFields): array([[ 2.]]) """ - mask = self.getMask(cutExpression, coordSys=coordSys) + mask = self.evalf(expression, coordSys=coordSys) inst = self[mask].copy() - inst.triangleScalars = inst.triangleScalars[mask.reshape((self.getNTriangles(), 3)).T[0]] + inst.triangleScalars = inst.triangleScalars[mask.reshape((self.get_ntriangles(), 3)).T[0]] return inst @decoTools.cached_property() - def areas(self): + def _areas(self): """ Cached method to retrieve areas of triangles """ - indices = range(self.getNTriangles()) - aIndices = [i * 3 for i in indices] - bIndices = [i * 3 + 1 for i in indices] - cIndices = [i * 3 + 2 for i in indices] - - # define 3 vectors building the triangle and take their norm - 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) - lengths.sort() - a, b, c = lengths.T + transform = np.eye() + return self.areas(transform=transform) - return 0.25 * np.sqrt((a + (b + c)) * (c - (a - b)) * (c + (a - b)) * (a + (b - c))) - - def getAreas(self): + def areas(self, transform=None): """ Calculate area with "heron's formula" + Args: + transform (np.ndarray): optional transformation matrix + The triangle points are transformed with transform if given + before calclulating the area Examples: >>> m = tfields.Mesh3D([[1,0,0], [0,0,1], [0,0,0]], [[0, 1, 2]]); - >>> m.triangles.getAreas() + >>> m.triangles.areas() array([ 0.5]) >>> m = tfields.Mesh3D([[1,0,0], [0,1,0], [0,0,0], [0,0,1]], [[0, 1, 2], [1, 2, 3]]); - >>> m.triangles.getAreas() + >>> m.triangles.areas() array([ 0.5, 0.5]) >>> m = tfields.Mesh3D([[1,0,0], [0,1,0], [1,1,0], [0,0,1], [1,0,1]], ... [[0, 1, 2], [0, 3, 4]]); - >>> m.triangles.getAreas() + >>> m.triangles.areas() array([ 0.5, 0.5]) """ - return self.areas - - def getAreasWithTransformation(self, transformationMatrix): - """ - Calculate area with "heron's formula". - The triangle points are transformed with transformationMatrix - """ - indices = range(self.getNTriangles()) - aIndices = [i * 3 for i in indices] - bIndices = [i * 3 + 1 for i in indices] - cIndices = [i * 3 + 2 for i in indices] - - # define 3 vectors building the triangle, transform it back and take their norm - - a = np.linalg.norm(np.linalg.solve(transformationMatrix.T, (self[aIndices, :] - self[bIndices, :]).T), axis=0) - b = np.linalg.norm(np.linalg.solve(transformationMatrix.T, (self[aIndices, :] - self[cIndices, :]).T), axis=0) - c = np.linalg.norm(np.linalg.solve(transformationMatrix.T, (self[bIndices, :] - self[cIndices, :]).T), axis=0) - - # sort by length for numerical stability - lengths = np.concatenate((a.reshape(-1, 1), b.reshape(-1, 1), c.reshape(-1, 1)), axis=1) - lengths.sort() - a, b, c = lengths.T - - return 0.25 * np.sqrt((a + (b + c)) * (c - (a - b)) * (c + (a - b)) * (a + (b - c))) - - def getCorners(self): + if transform is None: + return self._areas + else: + indices = range(self.get_ntriangles()) + aIndices = [i * 3 for i in indices] + bIndices = [i * 3 + 1 for i in indices] + cIndices = [i * 3 + 2 for i in indices] + + # define 3 vectors building the triangle, transform it back and take their norm + + if not np.array_equal(transform, np.eye): + 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) + + # sort by length for numerical stability + lengths = np.concatenate((a.reshape(-1, 1), b.reshape(-1, 1), c.reshape(-1, 1)), axis=1) + lengths.sort() + a, b, c = lengths.T + + return 0.25 * np.sqrt((a + (b + c)) * (c - (a - b)) * (c + (a - b)) * (a + (b - c))) + + def corners(self): """ Returns: three np.arrays with corner points of triangles """ - indices = range(self.getNTriangles()) + indices = range(self.get_ntriangles()) aIndices = [i * 3 for i in indices] bIndices = [i * 3 + 1 for i in indices] cIndices = [i * 3 + 2 for i in indices] @@ -215,21 +185,21 @@ class Triangles3D(tfields.TensorFields): c = self[cIndices, :] return a, b, c - def getCircumcenters(self): + def circumcenters(self): """ Semi baricentric method to calculate circumcenter points of the triangles Examples: >>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [-1,0,0], [0,1,0], [0,0,1]], ... [[0, 1, 3],[0, 2, 3],[1,2,4], [1, 3, 4]]); - >>> m.triangles.getCircumcenters() + >>> m.triangles.circumcenters() Points3D([[ 5.00000000e-01, 5.00000000e-01, 0.00000000e+00], [ -5.00000000e-01, 5.00000000e-01, 0.00000000e+00], [ 0.00000000e+00, 0.00000000e+00, 2.22044605e-16], [ 3.33333333e-01, 3.33333333e-01, 3.33333333e-01]]) """ - pointA, pointB, pointC = self.getCorners() + pointA, pointB, pointC = self.corners() 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) @@ -246,134 +216,115 @@ class Triangles3D(tfields.TensorFields): return tfields.Points3D(P) @decoTools.cached_property() - def centroids(self): + def _centroids(self): """ this version is faster but takes much more ram also. So much that i get memory error with a 32 GB RAM """ - nT = self.getNTriangles() + nT = self.get_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], coordSys=self.coordSys) """ Old version: - pointA, pointB, pointC = self.getCorners() + pointA, pointB, pointC = self.corners() return Points3D(1. / 3 * (pointA + pointB + pointC)), coordSys=self.coordSys This versioin was slightly slower (110 % of new version) Memory usage of new version is better for a factor of 4 or so. Not really reliable method of measurement """ - def getCentroids(self): + def centroids(self): """ Examples: >>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [-1,0,0], [0,1,0], [0,0,1]], ... [[0, 1, 3],[0, 2, 3],[1,2,4], [1, 3, 4]]); - >>> m.triangles.getCentroids() + >>> m.triangles.centroids() Points3D([[ 0.33333333, 0.33333333, 0. ], [-0.33333333, 0.33333333, 0. ], [ 0. , 0. , 0.33333333], [ 0.33333333, 0.33333333, 0.33333333]]) """ - return self.centroids + return self._centroids - def getTriangleCenters(self): - """ - Deprecated! use centroids + def edges(self): """ - warnings.warn("deprecated", DeprecationWarning) - return self.centroids - - def getEdgeVectors(self): - """ - Return two vectors that define two triangle edges + Retrieve two of the three edge vectors + Returns: + two np.ndarrays: vectors ab and ac, where a, b, c are corners (see + self.corners) """ - a, b, c = self.getCorners() + a, b, c = self.corners() ab = b - a ac = c - a return ab, ac - def getNormVectors(self): + def norms(self): """ Examples: >>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [-1,0,0], [0,1,0], [0,0,1]], ... [[0, 1, 3],[0, 2, 3],[1,2,4], [1, 3, 4]]); - >>> m.triangles.getNormVectors() + >>> m.triangles.norms() array([[ 0. , 0. , 1. ], [ 0. , 0. , -1. ], [ 0. , 1. , 0. ], [ 0.57735027, 0.57735027, 0.57735027]]) """ - ab, ac = self.getEdgeVectors() + ab, ac = self.edges() vectors = np.cross(ab, ac) norms = np.apply_along_axis(np.linalg.norm, 0, vectors.T).reshape(-1, 1) # cross product may be zero, so replace zero norms by ones to divide vectors by norms np.place(norms, norms == 0.0, 1.0) return vectors / norms - def getBaricentricCoordinates(self, point, delta=0.): + def _baricentric(self, point, delta=0.): """ Determine baricentric coordinates like [u,v,w] = [ab, ac, ab x ac]^-1 * ap where ax is vector from point a to point x - return baricentric coordinates of triangle - last coordinate is prism Examples: empty Meshes return right formatted array >>> m = tfields.Mesh3D([], faces=[]) - >>> m.triangles.getBaricentricCoordinates(np.array([0.2, 0.2, 0])) + >>> m.triangles._baricentric(np.array([0.2, 0.2, 0])) array([], dtype=float64) - >>> m2 = tfields.Mesh3D([[0,0,0], [2,0,0], [0,2,0], [0,0,2]], [[0, 1, 2], [0, 2, 3]]); - >>> m2.triangles.getBaricentricCoordinates(np.array([0.2, 0.2, 0]), delta=2.) - array([[ 0.1, 0.1, 0. ], - [ 0.1, 0. , 0.1]]) + >>> m2 = tfields.Mesh3D([[0,0,0], [2,0,0], [0,2,0], [0,0,2]], + ... faces=[[0, 1, 2], [0, 2, 3]]); + >>> assert np.array_equal(m2.triangles._baricentric(np.array([0.2, 0.2, 0]), + ... delta=2.), + ... [[0.1, 0.1, 0.0], + ... [0.1, 0.0, 0.1]]) - if point lies in the plane, return nan, else inf for w if delta is exactly 0. - >>> m2.triangles.getBaricentricCoordinates(np.array([0.2, 0.2, 0]), delta=0.) - array([[ 0.1, 0.1, nan], - [ 0.1, 0. , inf]]) + if point lies in the plane, return np.nan, else inf for w if delta is exactly 0. + >>> assert np.array_equal(m2.triangles._baricentric(np.array([0.2, 0.2, 0]), + ... delta=0.), + ... [[0.1, 0.1, np.nan], + ... [0.1, 0.0, np.inf]]) If you define triangles that have colinear side vectors or in general lead to - not invertable matrices [ab, ac, ab x ac] a warning wil be thrown but no error + not invertable matrices [ab, ac, ab x ac] the values will be nan >>> m3 = tfields.Mesh3D([[0,0,0], [2,0,0], [4,0,0], [0,1,0]], [[0, 1, 2], [0, 1, 3]]); - >>> import pytest - >>> with pytest.warns(UserWarning) as record: - ... bc = m3.triangles.getBaricentricCoordinates(np.array([0.2, 0.2, 0]), delta=0.3) - >>> bc - array([[ nan, nan, nan], - [ 0.1, 0.2, 0. ]]) - >>> len(record) - 1 - >>> "Could not invert matrix" in str(record[0].message) - True + >>> bc = m3.triangles._baricentric(np.array([0.2, 0.2, 0]), delta=0.3) + >>> assert np.array_equal(bc, [[ np.nan, np.nan, np.nan], [ 0.1, 0.2, 0. ]]) Returns: - arrays of u, v, w for each triangle + np.ndarray: barycentric coordinates u, v, w of point with respect to each triangle """ - if self.getNTriangles() == 0: + if self.get_ntriangles() == 0: return np.array([]) - 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)) - except: - raise - - a, b, c = self.getCorners() + a, _, _ = self.corners() ap = point - a # matrix vector product for matrices and vectors - barCoords = np.einsum('...ji,...i', self.barycentricMatrix, ap) + barCoords = np.einsum('...ji,...i', self._baricentric_matrix, ap) with warnings.catch_warnings(): warnings.filterwarnings('ignore', message="invalid value encountered in divide") warnings.filterwarnings('ignore', message="divide by zero encountered in divide") @@ -381,25 +332,22 @@ class Triangles3D(tfields.TensorFields): return barCoords @decoTools.cached_property() - def barycentricMatrix(self): + def _baricentric_matrix(self): """ cached barycentric matrix for faster calculations """ - a, b, c = self.getCorners() - ab = b - a - ac = c - a + ab, ac = self.edges() - # get norm vector + # get norm vector TODO: replace by norm = self.norms() norm = np.cross(ab, ac) normLen = np.linalg.norm(norm, axis=1) normLen = normLen.reshape((1,) + normLen.shape) - colinearMask = (normLen == 0) + colinear_mask = (normLen == 0) with warnings.catch_warnings(): warnings.filterwarnings('ignore', category=np.VisibleDeprecationWarning) # prevent divide by 0 - colinSlice = np.where(~colinearMask.T) - norm[np.where(~colinearMask.T)] = \ - norm[np.where(~colinearMask.T)] / normLen.T[np.where(~colinearMask.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 @@ -418,108 +366,21 @@ class Triangles3D(tfields.TensorFields): matrixI.append(np.full((3, 3), np.nan)) return np.array(matrixI) - def pointsInTriangles(self, points3D, delta=0., method='baryc', assignMultiple=False): + def _in_triangles(self, point, delta=0.): """ - Args: - points3D (Points3D instance) - - optional: - delta (float): normal distance to a triangle, that the points - is concidered to be contained in the triangle. - method (str): choose, which method to use for the calculation. - baryc is recommended - assignMultiple (bool): 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: - 2-d numpy mask which is True, where a point is in a face - face ind -> - points index - 1 0 0 - | 1 0 0 - V 0 0 1 - - if assignMultiple == False: - there is just one True for each points index - face indices can have multiple true values - - """ - log = logger.new() - log.verbose("Checking for points in Triangles") - try: - points3D = points3D.reshape((-1, 3)) - except ValueError: - raise ValueError("points3D must be castable to shape (nPoints, 3) " - "but is of shape {0}".format(points3D.shape)) - except: - raise - - if self.getNTriangles() == 0: - return np.empty((points3D.shape[0], 0), dtype=bool) - - masks = np.zeros((points3D.getNPoints(), self.getNTriangles()), dtype=bool) - for i, point in enumerate(points3D): - if i % 1000 == 0: - log.verbose("Status: {i} points processed. Yet {nP} are found " - "to be in triangles.".format(nP=masks.sum() - 1, **locals())) - masks[i] = self.pointInTriangles(point, delta, method=method) - log.verbose("Status: All points processed. Yet {nP} are found " - "to be in triangles.".format(nP=masks.sum() - 1, **locals())) - - if len(masks) == 0: - return masks - - if not assignMultiple: - """ - 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 a matrix: - face ind -> - points index - 1 0 0 - | 1 0 0 - V 0 0 1 - - there is just one True for each points index - face indices can have multiple true values - """ - return masks - - def pointInTriangles(self, point, delta=0., method='baryc'): - """ - Returns: - np.array: boolean mask, True where point in a triangle within delta + Barycentric method to optain, wheter a point is in any of the triangles Args: point (list of len 3) delta (float): acceptance in +- norm vector direction - """ - if method == 'baryc': - return self.pointInTrianglesBaryc(point, delta=delta) - elif method == 'baryc2': - return self.pointInTrianglesBarycentric2(point, delta=delta) - else: - raise NotImplementedError - - def pointInTrianglesBaryc(self, point, delta=0.): - """ + Returns: + 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]], [[0, 1, 2]]); - >>> m.triangles.pointInTrianglesBaryc(tfields.Points3D([[0.2, 0.2, 0]])) + >>> m.triangles._in_triangles(tfields.Points3D([[0.2, 0.2, 0]])) array([ True], dtype=bool) wrong format of point will raise a ValueError - >>> m.triangles.pointInTrianglesBaryc(tfields.Points3D([[0.2, 0.2, 0], [0.2, 0.2, 0]])) + >>> 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) @@ -528,15 +389,15 @@ class Triangles3D(tfields.TensorFields): >>> m2 = tfields.Mesh3D([[1,0,0], [0,1,0], [0,0,0], [4,0,0], [4, 4, 0], [8, 0, 0]], ... [[0, 1, 2], [3, 4, 5]]); - >>> m2.triangles.pointInTrianglesBaryc(np.array([0.2, 0.2, 0])) + >>> m2.triangles._in_triangles(np.array([0.2, 0.2, 0])) array([ True, False], dtype=bool) - >>> m2.triangles.pointInTrianglesBaryc(np.array([5, 2, 0])) + >>> m2.triangles._in_triangles(np.array([5, 2, 0])) array([False, True], dtype=bool) delta allows to accept points that lie within delta orthogonal to the tringle plain - >>> m2.triangles.pointInTrianglesBaryc(np.array([0.2, 0.2, 9000]), 0.0) + >>> m2.triangles._in_triangles(np.array([0.2, 0.2, 9000]), 0.0) array([False, False], dtype=bool) - >>> m2.triangles.pointInTrianglesBaryc(np.array([0.2, 0.2, 0.1]), 0.2) + >>> m2.triangles._in_triangles(np.array([0.2, 0.2, 0.1]), 0.2) array([ True, False], dtype=bool) If you define triangles that have colinear side vectors or in general lead to @@ -544,93 +405,91 @@ class Triangles3D(tfields.TensorFields): >>> m3 = tfields.Mesh3D([[0,0,0], [2,0,0], [4,0,0], [0,1,0]], [[0, 1, 2], [0, 1, 3]]); >>> import pytest >>> with pytest.warns(UserWarning) as record: - ... mask = m3.triangles.pointInTrianglesBaryc(np.array([0.2, 0.2, 0]), delta=0.3) + ... mask = m3.triangles._in_triangles(np.array([0.2, 0.2, 0]), delta=0.3) >>> mask array([False, True], dtype=bool) """ - if self.getNTriangles() == 0: + if self.get_ntriangles() == 0: return np.array([], dtype=bool) - u, v, w = self.getBaricentricCoordinates(point, delta=delta).T + u, v, w = self._baricentric(point, delta=delta).T if delta == 0.: - w[np.isnan(w)] = 0. # division by 0 in getBaricentricCoordinates makes w = 0 nan. + w[np.isnan(w)] = 0. # division by 0 in baricentric makes w = 0 nan. with warnings.catch_warnings(): warnings.filterwarnings('ignore', message="invalid value encountered in less_equal") orthogonalAcceptance = (abs(w) <= 1) - barycentricBools = ((v <= 1.) & (v >= 0.)) & ((u <= 1.) & (u >= 0.)) & ((v + u <= 1.0)) - return np.array(barycentricBools & orthogonalAcceptance) + barycentric_bools = ((v <= 1.) & (v >= 0.)) & ((u <= 1.) & (u >= 0.)) & ((v + u <= 1.0)) + return np.array(barycentric_bools & orthogonalAcceptance) - def pointInTrianglesBarycentric2(self, point, delta=0.): + def in_triangles(self, tensors, delta=0., assign_multiple=False): """ - Check whether point lies within triangles with Barycentric Technique - It could lie within multiple triangles! - Examples: - >>> m = tfields.Mesh3D([[1,0,0], [0,1,0], [0,0,0]], [[0, 1, 2]]); - >>> m.triangles.pointInTrianglesBarycentric2(tfields.Points3D([[0.2, 0.2, 0]])) - array([ True], dtype=bool) + Barycentric method to obtain, which tensors are containes in any of the triangles + Args: + tensors (Points3D instance) - wrong format of point will raise a ValueError - >>> m.triangles.pointInTrianglesBarycentric2(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) + optional: + delta (float): normal distance to a triangle, that the points + is concidered to be contained in the triangle. + assign_multiple (bool): 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 - 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]], - ... [[0, 1, 2], [3, 4, 5]]); - >>> m2.triangles.pointInTrianglesBarycentric2(np.array([0.2, 0.2, 0])) - array([ True, False], dtype=bool) - >>> m2.triangles.pointInTrianglesBarycentric2(np.array([5, 2, 0])) - array([False, True], dtype=bool) + Returns: + np.ndarray: 2-d mask which is True, where a point is in a triangle + triangle index -> + points index + 1 0 0 + | 1 0 0 + V 0 0 1 - delta allows to accept points that lie within delta orthogonal to the tringle plain - >>> m2.triangles.pointInTrianglesBarycentric2(np.array([0.2, 0.2, 9000]), 0.0) - array([False, False], dtype=bool) - >>> m2.triangles.pointInTrianglesBarycentric2(np.array([0.2, 0.2, 0.1]), 0.2) - array([ True, False], dtype=bool) + if assign_multiple == False: + there is just one True for each points index + triangle indices can have multiple true values - Returns: - array of bools for each triangle """ - 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)) - except: - raise - - indices = range(self.getNTriangles()) - aIndices = [i * 3 for i in indices] - bIndices = [i * 3 + 1 for i in indices] - cIndices = [i * 3 + 2 for i in indices] + log = logging.getLogger(__name__) + if self.get_ntriangles() == 0: + return np.empty((tensors.shape[0], 0), dtype=bool) - # compute base vectors - AP = (self[aIndices, :] - point).T - AB = (self[aIndices, :] - self[bIndices, :]).T - AC = (self[aIndices, :] - self[cIndices, :]).T + masks = np.zeros((len(tensors), self.get_ntriangles()), dtype=bool) + with tensors.tmp_transform(self.coordSys): + for i, point in enumerate(tensors): + if i % 1000 == 0: + log.debug("Status: {i} points processed. Yet {nP} are found " + "to be in triangles.".format(nP=masks.sum() - 1, **locals())) + masks[i] = self._in_triangles(point, delta) + log.debug("Status: All points processed. Yet {nP} are found " + "to be in triangles.".format(nP=masks.sum() - 1, **locals())) - # compute barycentric coordinates - np.seterr(divide='ignore', invalid='ignore') + if not masks: + # empty sequence + return masks - denom = AB[1] * AC[0] - AB[0] * AC[1] - v = (AB[1] * AP[0] - AB[0] * AP[1]) / denom - u = -(AC[1] * AP[0] - AC[0] * AP[1]) / denom - # u = (AP[0] - v * AC[0]) / AB[0] # also true, but problematic with AB[0] + 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) - # calculate bools - orthogonalAcceptance = abs(AP[2] - u * AB[2] - v * AC[2]) <= delta - barycentricBools = ((v <= 1.) & (v >= 0.)) & ((u <= 1.) & (u >= 0.)) & ((v + u <= 1.0)) - np.seterr(divide='warn', invalid='warn') - return np.array(barycentricBools & orthogonalAcceptance) + 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 - def pointOnTriangleSide(self, point): + def _on_edges(self, point): """ + TODO: on_edges like in_triangles + Determine whether a point is on the edge / side ray of a triangle Returns: np.array: boolean mask which is true, if point is on one side ray of a triangle @@ -639,157 +498,40 @@ class Triangles3D(tfields.TensorFields): ... [[0, 1, 3],[0, 2, 3],[1,2,4]]); Corner points are found - >>> m.triangles.pointOnTriangleSide(tfields.Points3D([[0,1,0]])) + >>> m.triangles.on_edges(tfields.Points3D([[0,1,0]])) array([ True, True, False], dtype=bool) Side points are found, too - >>> m.triangles.pointOnTriangleSide(tfields.Points3D([[0.5,0,0.5]])) + >>> m.triangles.on_edges(tfields.Points3D([[0.5,0,0.5]])) array([False, False, True], dtype=bool) """ - u, v, w = self.getBaricentricCoordinates(point, 1.).T + u, v, w = self._baricentric(point, 1.).T orthogonalAcceptance = (w == 0) # point should lie in triangle - barycentricBools = (((0. <= v) & (v <= 1.)) & (u == 0.)) | \ + barycentric_bools = (((0. <= v) & (v <= 1.)) & (u == 0.)) | \ (((0. <= u) & (u <= 1.)) & (v == 0.)) | \ (v + u == 1.0) - return np.array(barycentricBools & orthogonalAcceptance) - - def getWeightedAreas(self, weightScalarIndex=0): - """ - Args: - weightScalarIndex (int): index to scalars to weight with - """ - # set weights to 1.0 if weightScalarIndex is None - if weightScalarIndex is None: - weights = np.ones(self.getNTriangles()) - elif self.getScalarArrays().shape[0] == 0: - weights = np.ones(self.getNTriangles()) - else: - weights = self.getScalarArrays()[weightScalarIndex] - areas = self.getAreas() - return areas * weights - - def getNormedWeightedAreas(self, weightScalarIndex=0, transform=None): - """ - Args: - weightScalarIndex (int): index to scalars to weight with - Returns: - A_i * w_i / sum_{i=0}^{N}(A_i * w_i) - """ - # set weights to 1.0 if weightScalarIndex is None - if transform is None: - transform = np.eye(3) - if weightScalarIndex is None: - weights = np.ones(self.getNTriangles()) - elif self.getScalarArrays().shape[0] == 0: - weights = np.ones(self.getNTriangles()) - else: - weights = self.getScalarArrays()[weightScalarIndex] - areas = self.getAreasWithTransformation(transform) - return areas * weights / np.dot(areas, weights) + return np.array(barycentric_bools & orthogonalAcceptance) - def getMean(self, weightScalarIndex=0): + def _weights(self, weights, rigid=False): """ + transformer method for weights inputs. Args: - weightScalarIndex (int): index to scalars to weight with + weights (np.ndarray | int | None): + If weights is integer it will be used as index for fields and fields are + used as weights. + If weights is None it will + Otherwise just pass the weights. Returns: - mean of the centroids weighted by scalars @ <weightScalarIndex> - """ - centroids = self.getCentroids() - weightedAreas = self.getNormedWeightedAreas(weightScalarIndex=weightScalarIndex) - - # sum up all triangle centers weighted by the area and divided by the total area - return np.dot(centroids.T, weightedAreas) - def getStd(self, mean=None, weightScalarIndex=0): """ - Args: - mean(Points3D) - weightScalarIndex (int): index to scalars to weight with - """ - if mean is None: - mean = self.getMean(weightScalarIndex=weightScalarIndex) - centroids = self.getCentroids() - weightedAreas = self.getNormedWeightedAreas(weightScalarIndex=weightScalarIndex) - squaredDiff = (centroids - mean) ** 2 - - var = np.dot(squaredDiff.T, weightedAreas) - return np.sqrt(var) - - def getDirection(self, weightScalarIndex=0, transform=None): - """ - Returns: Eigenvector corresponding to greatest eigenvalue multiplied by greatest eigenvalue - """ - # calculate covariance matrix - log = logger.new() - try: - cov = np.cov(self.getCentroids().T, - ddof=0, - aweights=self.getNormedWeightedAreas(weightScalarIndex=weightScalarIndex, transform=transform)) - # calculate eigenvalues and eigenvectors of covariance - - evals, evecs = np.linalg.eigh(cov) - except: - log.warning("np.linalg.eig() failed. use zeros instead!") - evals = np.zeros(3) - evecs = np.zeros((3, 3)) - idx = np.argmax(evals) - # return eigenvector corresponding to greatest eigenvalue multiplied by greatest eigenvalue - # always positive in x-direction - return evals.max() * (evecs[:, idx] if evecs[0, idx] > 0 else -evecs[:, idx]) - - def getCovEig(self, weightScalarIndex=0): - """ - Args: - weightScalarIndex (int): index to scalars to weight with - """ - cov = np.cov(self.getCentroids().T, - ddof=0, - aweights=self.getNormedWeightedAreas(weightScalarIndex=weightScalarIndex)) - # calculate eigenvalues and eigenvectors of covariance - evals, evecs = np.linalg.eigh(cov) - idx = evals.argsort()[::-1] - evals = evals[idx] - evecs = evecs[:, idx] - e = np.concatenate((evecs, evals.reshape(1, 3))) - return e.T.reshape(12, ) - - def getMainAxes(self): - """ - Args: - weightScalarIndex (int): index to scalars to weight with - Returns: - Main Axes eigen-vectors - """ - mean = self.getMean(None) - relativeCentroids = self.getCentroids() - mean - x = relativeCentroids - cov = np.cov(x.T, - ddof=0, - aweights=self.getNormedWeightedAreas(weightScalarIndex=None)) - # calculate eigenvalues and eigenvectors of covariance - evals, evecs = np.linalg.eigh(cov) - return (evecs * evals.T).T - - def getMesh3D(self): - """ - Forward to an alternative representation of the triangeles as mesh - """ - return tfields.Mesh3D(self, - [[3 * i, 3 * i + 1, 3 * i + 2] - for i in range(self.getNTriangles())], - faceScalars=self.triangleScalars) - - def plot(self, **kwargs): - """ - Forward to Mesh3D plot - """ - return self.getMesh3D().plot(**kwargs) + # set weights to 1.0 if weights is None + if weights is None: + weights = np.ones(self.get_ntriangles()) + return super(Triangles3D, self)._weights(weights, rigid=rigid) if __name__ == '__main__': import doctest - import sympy # NOQA: F401 - doctest.testmod() -- GitLab