From 08eb56916bd0357a24b599c098c3e868a21b8a95 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Daniel=20B=C3=B6ckenhoff=20=28Laptop=29?= <dboe@ipp.mpg.de>
Date: Sat, 19 May 2018 16:33:07 +0200
Subject: [PATCH] added array_Set_ops, bases as file
---
tfields/array_set_ops.py | 107 ++++++
tfields/bases/__init__.py | 151 +++++++++
tfields/bases/manifold_3.py | 83 +++++
tfields/points3D.py | 630 ++++++++++++++++++++++++++++++++++++
4 files changed, 971 insertions(+)
create mode 100644 tfields/array_set_ops.py
create mode 100644 tfields/bases/__init__.py
create mode 100644 tfields/bases/manifold_3.py
create mode 100644 tfields/points3D.py
diff --git a/tfields/array_set_ops.py b/tfields/array_set_ops.py
new file mode 100644
index 0000000..54025ac
--- /dev/null
+++ b/tfields/array_set_ops.py
@@ -0,0 +1,107 @@
+#!/usr/bin/env
+# encoding: utf-8
+"""
+Author: Daniel Boeckenhoff
+Mail: daniel.boeckenhoff@ipp.mpg.de
+
+Collection of additional numpy functions
+part of tfields library
+"""
+import numpy as np
+import loggingTools
+logger = loggingTools.Logger(__name__)
+
+
+def convert_nan(array, value=0.):
+ """
+ Replace all occuring NaN values by value
+ """
+ nanIndices = np.isnan(array)
+ array[nanIndices] = value
+
+
+def view1D(ar):
+ """
+ https://stackoverflow.com/a/44999009/ @Divakar
+ """
+ ar = np.ascontiguousarray(ar)
+ voidDt = np.dtype((np.void, ar.dtype.itemsize * ar.shape[1]))
+ return ar.view(voidDt).ravel()
+
+
+def argsort_unique(idx):
+ """
+ https://stackoverflow.com/a/43411559/ @Divakar
+ """
+ n = idx.size
+ sidx = np.empty(n, dtype=int)
+ sidx[idx] = np.arange(n)
+ return sidx
+
+
+def duplicates(ar, axis=None):
+ """
+ View1D version of duplicate search
+ Speed up version after
+ https://stackoverflow.com/questions/46284660 \
+ /python-numpy-speed-up-2d-duplicate-search/46294916#46294916
+ Args:
+ ar (array_like): array
+ other args: see np.isclose
+ Examples:
+ >>> import tfields
+ >>> import numpy as np
+ >>> a = np.array([[1, 0, 0], [1, 0, 0], [2, 3, 4]])
+ >>> tfields.duplicates(a, axis=0)
+ array([0, 0, 2])
+
+ Returns:
+ list of int: int is pointing to first occurence of unique value
+ """
+ if axis != 0:
+ raise NotImplementedError()
+ sidx = np.lexsort(ar.T)
+ b = ar[sidx]
+
+ groupIndex0 = np.flatnonzero((b[1:] != b[:-1]).any(1)) + 1
+ groupIndex = np.concatenate(([0], groupIndex0, [b.shape[0]]))
+ ids = np.repeat(range(len(groupIndex) - 1), np.diff(groupIndex))
+ sidx_mapped = argsort_unique(sidx)
+ ids_mapped = ids[sidx_mapped]
+
+ grp_minidx = sidx[groupIndex[:-1]]
+ out = grp_minidx[ids_mapped]
+ return out
+
+
+def index(ar, entry, rtol=0, atol=0, equal_nan=False, axis=None):
+ """
+ Examples:
+ >>> import tfields
+ >>> a = np.array([[1, 0, 0], [1, 0, 0], [2, 3, 4]])
+ >>> tfields.index(a, [2, 3, 4], axis=0)
+ 2
+
+ >>> a = np.array([[1, 0, 0], [2, 3, 4]])
+ >>> tfields.index(a, 4)
+ 5
+
+ Returns:
+ list of int: indices of point occuring
+ """
+ if axis is None:
+ ar = ar.flatten()
+ elif axis != 0:
+ raise NotImplementedError()
+ for i, part in enumerate(ar):
+ isclose = np.isclose(part, entry, rtol=rtol, atol=atol,
+ equal_nan=equal_nan)
+ if axis is not None:
+ isclose = isclose.all()
+ if isclose:
+ return i
+
+
+if __name__ == '__main__':
+ import doctest
+ doctest.testmod()
diff --git a/tfields/bases/__init__.py b/tfields/bases/__init__.py
new file mode 100644
index 0000000..e45d03e
--- /dev/null
+++ b/tfields/bases/__init__.py
@@ -0,0 +1,151 @@
+#!/usr/bin/env
+# encoding: utf-8
+"""
+Author: Daniel Boeckenhoff
+Mail: daniel.boeckenhoff@ipp.mpg.de
+
+part of tfields library
+Tools for sympy coordinate transformation
+"""
+import tfields
+import numpy as np
+import sympy
+import sympy.diffgeom
+from six import string_types
+import warnings
+
+from .manifold_3 import CARTESIAN, CYLINDER, SPHERICAL
+from .manifold_3 import cartesian, cylinder, spherical
+
+
+def get_coord_system(base):
+ """
+ Args:
+ base (str or sympy.diffgeom.get_coord_system)
+ Return:
+ sympy.diffgeom.get_coord_system
+ """
+ if isinstance(base, string_types):
+ base = getattr(tfields.bases, base)
+ if not isinstance(base, sympy.diffgeom.CoordSystem):
+ raise TypeError("Wrong type of coordSystem base.")
+ return base
+
+
+def get_coord_system_name(base):
+ """
+ Args:
+ base (str or sympy.diffgeom.get_coord_system)
+ Returns:
+ str: name of base
+ """
+ if isinstance(base, sympy.diffgeom.CoordSystem):
+ base = str(getattr(base, 'name'))
+ # if not (isinstance(base, string_types) or base is None):
+ # baseType = type(base)
+ # raise ValueError("Coordinate system must be string_type."
+ # " Retrieved value '{base}' of type {baseType}."
+ # .format(**locals()))
+ return base
+
+
+def lambdifiedTrafo(base_old, base_new):
+ """
+ Args:
+ base_old (sympy.CoordSystem)
+ base_new (sympy.CoordSystem)
+
+ Examples:
+ >>> import numpy as np
+ >>> import tfields
+
+ Transform cartestian to cylinder or spherical
+ >>> a = np.array([[3,4,0]])
+
+ >>> trafo = tfields.bases.lambdifiedTrafo(tfields.bases.cartesian,
+ ... tfields.bases.cylinder)
+ >>> new = np.concatenate([trafo(*coords).T for coords in a])
+ >>> assert new[0, 0] == 5
+
+ >>> trafo = tfields.bases.lambdifiedTrafo(tfields.bases.cartesian,
+ ... tfields.bases.spherical)
+ >>> new = np.concatenate([trafo(*coords).T for coords in a])
+ >>> assert new[0, 0] == 5
+
+ """
+ coords = tuple(base_old.coord_function(i) for i in range(base_old.dim))
+ f = sympy.lambdify(coords,
+ base_old.coord_tuple_transform_to(base_new,
+ list(coords)),
+ modules='numpy')
+ return f
+
+
+def transform(array, base_old, base_new):
+ """
+ Transform the input array in place
+ Args:
+ array (np.ndarray)
+ base_old (str or sympy.CoordSystem):
+ base_new (str or sympy.CoordSystem):
+ Examples:
+ Cylindrical coordinates
+ >>> import tfields
+ >>> print "TEST"
+ >>> cart = np.array([[0, 0, 0],
+ ... [1, 0, 0],
+ ... [1, 1, 0],
+ ... [0, 1, 0],
+ ... [-1, 1, 0],
+ ... [-1, 0, 0],
+ ... [-1, -1, 0],
+ ... [0, -1, 0],
+ ... [1, -1, 0],
+ ... [0, 0, 1]])
+ >>> cyl = tfields.bases.transform(cart, 'cartesian', 'cylinder')
+ >>> cyl
+ >>> import npTools
+ >>> pnp = npTools.Points3D(cart)
+ >>> pnp.transform('cylinder')
+ >>> pnp
+
+ Transform cylinder to spherical. No connection is defined so routing via
+ cartesian
+ >>> import numpy as np
+ >>> import tfields
+ >>> b = np.array([[5, np.arctan(4. / 3), 0]])
+ >>> newB = tfields.bases.transform(b, 'cylinder', 'spherical')
+ >>> assert newB[0, 0] == 5
+ >>> assert round(newB[0, 1], 10) == round(b[0, 1], 10)
+
+ """
+ base_old = get_coord_system(base_old)
+ base_new = get_coord_system(base_new)
+ if base_new not in base_old.transforms:
+ for baseTmp in base_new.transforms:
+ if baseTmp in base_old.transforms:
+ transform(array, base_old, baseTmp)
+ transform(array, baseTmp, base_new)
+ return
+ raise ValueError("Trafo not found.")
+
+ # very fast trafos in numpy only
+ shortTrafo = None
+ try:
+ shortTrafo = getattr(base_old, 'to_{base_new.name}'.format(**locals()))
+ except AttributeError:
+ pass
+ if shortTrafo:
+ shortTrafo(array)
+ return
+
+ # trafo via lambdified sympy expressions
+ trafo = tfields.bases.lambdifiedTrafo(base_old, base_new)
+ with warnings.catch_warnings():
+ warnings.filterwarnings('ignore', message="invalid value encountered in double_scalars")
+ array[:] = np.concatenate([trafo(*coords).T for coords in array])
+
+
+if __name__ == '__main__': # pragma: no cover
+ import doctest
+ doctest.testmod()
diff --git a/tfields/bases/manifold_3.py b/tfields/bases/manifold_3.py
new file mode 100644
index 0000000..337f1fd
--- /dev/null
+++ b/tfields/bases/manifold_3.py
@@ -0,0 +1,83 @@
+import tfields
+import sympy
+import numpy as np
+
+CARTESIAN = 'cartesian'
+CYLINDER = 'cylinder'
+SPHERICAL = 'spherical'
+
+m = sympy.diffgeom.Manifold('M', 3)
+patch = sympy.diffgeom.Patch('P', m)
+
+# cartesian
+x, y, z = sympy.symbols('x, y, z')
+cartesian = sympy.diffgeom.CoordSystem(CARTESIAN, patch, ['x', 'y', 'z'])
+
+# cylinder
+
+R, Phi, Z = sympy.symbols('R, Phi, Z')
+cylinder = sympy.diffgeom.CoordSystem(CYLINDER, patch, ['R', 'Phi', 'Z'])
+cylinder.connect_to(cartesian,
+ [R, Phi, Z],
+ [R * sympy.cos(Phi),
+ R * sympy.sin(Phi),
+ Z],
+ inverse=False,
+ fill_in_gaps=False)
+cartesian.connect_to(cylinder,
+ [x, y, z],
+ [sympy.sqrt(x**2 + y**2),
+ sympy.Piecewise(
+ (sympy.pi, ((x < 0) & sympy.Eq(y, 0))),
+ (0., sympy.Eq(sympy.sqrt(x**2 + y**2), 0)),
+ (sympy.sign(y) * sympy.acos(x / sympy.sqrt(x**2 + y**2)), True)),
+ z],
+ inverse=False,
+ fill_in_gaps=False)
+
+
+def cylinder_to_cartesian(array):
+ rPoints = np.copy(array[:, 0])
+ phiPoints = np.copy(array[:, 1])
+ array[:, 0] = rPoints * np.cos(phiPoints)
+ array[:, 1] = rPoints * np.sin(phiPoints)
+
+
+def cartesian_to_cylinder(array):
+ x_vals = np.copy(array[:, 0])
+ y_vals = np.copy(array[:, 1])
+ problemPhiIndices = np.where((x_vals < 0) & (y_vals == 0))
+ array[:, 0] = np.sqrt(x_vals**2 + y_vals**2) # r
+ np.seterr(divide='ignore', invalid='ignore')
+ # phi
+ array[:, 1] = np.sign(y_vals) * np.arccos(x_vals / array[:, 0])
+ np.seterr(divide='warn', invalid='warn')
+ array[:, 1][problemPhiIndices] = np.pi
+ tfields.convert_nan(array, 0.) # for cases like cartesian 0, 0, 1
+
+
+cylinder.to_cartesian = cylinder_to_cartesian
+cartesian.to_cylinder = cartesian_to_cylinder
+
+# spherical
+r, phi, theta = sympy.symbols('r, phi, theta')
+spherical = sympy.diffgeom.CoordSystem(SPHERICAL, patch, ['r', 'phi', 'theta'])
+spherical.connect_to(cartesian,
+ [r, phi, theta],
+ [r * sympy.cos(phi) * sympy.sin(theta),
+ r * sympy.sin(phi) * sympy.sin(theta),
+ r * sympy.cos(theta)],
+ inverse=False,
+ fill_in_gaps=False)
+cartesian.connect_to(spherical,
+ [x, y, z],
+ [sympy.sqrt(x**2 + y**2 + z**2),
+ sympy.Piecewise(
+ (0., (sympy.Eq(x, 0) & sympy.Eq(y, 0))),
+ (sympy.atan(y / x), True)),
+ sympy.Piecewise(
+ (0., sympy.Eq(x**2 + y**2 + z**2, 0)),
+ # (0., (sympy.Eq(x, 0) & sympy.Eq(y, 0) & sympy.Eq(z, 0))),
+ (sympy.acos(z / sympy.sqrt(x**2 + y**2 + z**2)), True))],
+ inverse=False,
+ fill_in_gaps=False)
diff --git a/tfields/points3D.py b/tfields/points3D.py
new file mode 100644
index 0000000..862947b
--- /dev/null
+++ b/tfields/points3D.py
@@ -0,0 +1,630 @@
+#!/usr/bin/env
+# encoding: utf-8
+"""
+Author: Daniel Boeckenhoff
+Mail: daniel.boeckenhoff@ipp.mpg.de
+
+basic threedimensional tensor
+"""
+import tfields
+import numpy as np
+import os
+import osTools
+import ioTools
+import pyTools
+import mplTools as mpt
+import sympy
+import scipy as sp
+import scipy.spatial # NOQA: F401
+import scipy.stats as stats
+import warnings
+import loggingTools
+logger = loggingTools.Logger(__name__)
+
+
+class Points3D(tfields.Tensors):
+ # pylint: disable=R0904
+ """
+ Points3D is a general class for 3D Point operations and storage.
+ Points are stored in np.arrays of shape (len, 3).
+ Thus the three coordinates of the Points stay close.
+
+ Args:
+ points3DInstance -> copy constructor
+ [points3DInstance1, points3DInstance2, ...] -> coordSys are correctly treated
+ list of coordinates (see examples)
+
+ Kwargs:
+ coordSys (str):
+ Use tfields.bases.CARTESIAN -> x, y, z
+ Use tfields.bases.CYLINDER -> r, phi, z
+ Use tfields.bases.SPHERICAL -> r, phi, theta
+
+ Examples:
+ Initializing with 3 vectors
+ >>> import tfields
+ >>> import numpy as np
+ >>> p1 = tfields.Points3D([[1., 2., 3.], [4., 5., 6.], [1, 2, -6]])
+ >>> p1
+ Points3D([[ 1., 2., 3.],
+ [ 4., 5., 6.],
+ [ 1., 2., -6.]])
+
+ Initializing with list of coordinates
+ >>> p2 = tfields.Points3D(np.array([[1., 2., 3., 4, 5,],
+ ... [4., 5., 6., 7, 8],
+ ... [1, 2, -6, -1, 0]]).T)
+ >>> p2
+ Points3D([[ 1., 4., 1.],
+ [ 2., 5., 2.],
+ [ 3., 6., -6.],
+ [ 4., 7., -1.],
+ [ 5., 8., 0.]])
+ >>> p2.transform(tfields.bases.CYLINDER)
+ >>> p2
+ Points3D([[ 4.12310563, 1.32581766, 1. ],
+ [ 5.38516481, 1.19028995, 2. ],
+ [ 6.70820393, 1.10714872, -6. ],
+ [ 8.06225775, 1.05165021, -1. ],
+ [ 9.43398113, 1.01219701, 0. ]])
+
+ Copy constructor with one instance preserves coordSys of instance
+ >>> assert tfields.Points3D(p2).coordSys == p2.coordSys
+
+ Unless you specify other:
+ >>> tfields.Points3D(p2, coordSys=tfields.bases.CARTESIAN)
+ Points3D([[ 1., 4., 1.],
+ [ 2., 5., 2.],
+ [ 3., 6., -6.],
+ [ 4., 7., -1.],
+ [ 5., 8., 0.]])
+
+ Copy constructor with many instances chooses majority of coordinates
+ systems to avoid much transformation
+ >>> tfields.Points3D.merged(p1, p2, p1)
+ Points3D([[ 1., 2., 3.],
+ [ 4., 5., 6.],
+ [ 1., 2., -6.],
+ [ 1., 4., 1.],
+ [ 2., 5., 2.],
+ [ 3., 6., -6.],
+ [ 4., 7., -1.],
+ [ 5., 8., 0.],
+ [ 1., 2., 3.],
+ [ 4., 5., 6.],
+ [ 1., 2., -6.]])
+ >>> p1.transform(tfields.bases.CYLINDER)
+
+ ... unless specified other. Here it is specified
+ >>> tfields.Points3D.merged(p1, p2, coordSys=tfields.bases.CYLINDER)
+ Points3D([[ 2.23606798, 1.10714872, 3. ],
+ [ 6.40312424, 0.89605538, 6. ],
+ [ 2.23606798, 1.10714872, -6. ],
+ [ 4.12310563, 1.32581766, 1. ],
+ [ 5.38516481, 1.19028995, 2. ],
+ [ 6.70820393, 1.10714872, -6. ],
+ [ 8.06225775, 1.05165021, -1. ],
+ [ 9.43398113, 1.01219701, 0. ]])
+
+ Shape is always (..., 3)
+ >>> p = tfields.Points3D([[1., 2., 3.], [4., 5., 6.],
+ ... [1, 2, -6], [-5, -5, -5], [1,0,-1], [0,1,-1]])
+ >>> p.shape
+ (6, 3)
+
+ Empty array will create an ndarray of the form (0, 3)
+ >>> tfields.Points3D([])
+ Points3D([], shape=(0, 3), dtype=float64)
+
+ Use it as np.ndarrays -> masking etc. is inherited
+ >>> mask = np.array([True, False, True, False, False, True])
+ >>> mp = p[mask].copy()
+
+ Copy constructor
+ >>> mp
+ Points3D([[ 1., 2., 3.],
+ [ 1., 2., -6.],
+ [ 0., 1., -1.]])
+ >>> tfields.Points3D(mp)
+ Points3D([[ 1., 2., 3.],
+ [ 1., 2., -6.],
+ [ 0., 1., -1.]])
+
+ Coordinate system is implemented. Default is cartesian
+ >>> p.transform(tfields.bases.CYLINDER); p
+ Points3D([[ 2.23606798, 1.10714872, 3. ],
+ [ 6.40312424, 0.89605538, 6. ],
+ [ 2.23606798, 1.10714872, -6. ],
+ [ 7.07106781, -2.35619449, -5. ],
+ [ 1. , 0. , -1. ],
+ [ 1. , 1.57079633, -1. ]])
+ >>> p.transform(tfields.bases.CARTESIAN); p
+ Points3D([[ 1.000000e+00, 2.000000e+00, 3.000000e+00],
+ [ 4.000000e+00, 5.000000e+00, 6.000000e+00],
+ [ 1.000000e+00, 2.000000e+00, -6.000000e+00],
+ [-5.000000e+00, -5.000000e+00, -5.000000e+00],
+ [ 1.000000e+00, 0.000000e+00, -1.000000e+00],
+ [ 6.123234e-17, 1.000000e+00, -1.000000e+00]])
+ >>> p.getPhi()*180/3.1415
+ Points3D([ 63.43681974, 51.34170594, 63.43681974, -135.00398161,
+ 0. , 90.00265441])
+
+ Appending to the array with numpy.append
+ >>> a = tfields.Points3D([[1,2,3]])
+ >>> b = tfields.Points3D([[4,5,6]])
+ >>> c = np.append(a, b, 0)
+ >>> c
+ array([[1., 2., 3.],
+ [4., 5., 6.]])
+
+ """
+ def __new__(cls, tensors, **kwargs):
+ if not issubclass(type(tensors), Points3D):
+ kwargs['dim'] = 3
+ return super(Points3D, cls).__new__(cls, tensors, **kwargs)
+
+ def contains(self, other):
+ """
+ Inspired by a speed argument @
+ stackoverflow.com/questions/14766194/testing-whether-a-numpy-array-contains-a-given-row
+ Examples:
+ >>> p = tfields.Points3D([[1,2,3], [4,5,6], [6,7,8]])
+ >>> p.contains([4,5,6])
+ True
+
+ """
+ return any(np.equal(self, other).all(1))
+
+ def save(self, filePath, extension='npz', **kwargs):
+ """
+ save a tensors object.
+ filePath (str or buffer)
+ [extension (str): only needed if filePath is buffer]
+ """
+ if not osTools.isBuffer(filePath):
+ extension = ioTools.getExtension(filePath)
+ if extension not in ['npz', 'obj', 'txt', 'ply']:
+ raise NotImplementedError("Extension {0} not implemented.".format(extension))
+
+ # create from buffer, if file is tar file
+ if osTools.isTarPath(filePath):
+ raise NotImplementedError("Saving to tar file not intended "
+ "since we want to make the tar-writing outside.")
+
+ if extension == 'npz':
+ self.savez(filePath)
+ if extension == 'obj':
+ return self.saveObj(filePath, **kwargs)
+ if extension == 'txt':
+ self.saveTxt(filePath)
+ if extension == 'ply':
+ self.saveTxt(filePath)
+ else:
+ raise NotImplementedError("Do not understand the format of file {0}".format(filePath))
+
+ def savePly(self, filePath):
+ from plyfile import PlyData, PlyElement
+ with self.tmp_transform(tfields.bases.CARTESIAN):
+ vertices = np.array([tuple(x) for x in self],
+ dtype=[('x', '<f4'), ('y', '<f4'), ('z', '<f4')])
+ element = PlyElement.describe(vertices, 'vertex')
+
+ filePath = osTools.resolve(filePath)
+ PlyData([element]).write(filePath)
+
+ def saveTxt(self, filePath):
+ """
+ Save obj as .txt file
+ """
+ if self.coordSys != tfields.bases.CARTESIAN:
+ cpy = self.copy()
+ cpy.transform(tfields.bases.CARTESIAN)
+ else:
+ cpy = self
+ with ioTools.TextFile(filePath, 'w') as f:
+ f.writeMatrix(self, seperator=' ', lineBreak='\n')
+
+ def saveObj(self, filePath, **kwargs):
+ """
+ Save obj as wavefront/.obj file
+ """
+ obj = kwargs.pop('object', None)
+ group = kwargs.pop('group', None)
+
+ filePath = filePath.replace('.obj', '')
+ fileDir, fileName = os.path.split(filePath)
+ with open(osTools.resolve(filePath + '.obj'), 'w') as f:
+ f.write("# File saved with tensors Points3D.saveObj method\n\n")
+ if obj is not None:
+ f.write("o {0}\n".format(obj))
+ if group is not None:
+ f.write("g {0}\n".format(group))
+ for vertex in self:
+ f.write("v {v[0]} {v[1]} {v[2]}\n".format(v=vertex))
+
+ def savez(self, filePath):
+ """
+ Args:
+ filePath (open file or str/unicode): destination to save file to.
+ Examples:
+ >>> from tempfile import NamedTemporaryFile
+ >>> outFile = NamedTemporaryFile(suffix='.npz')
+ >>> p = tfields.Points3D([[1., 2., 3.], [4., 5., 6.], [1, 2, -6]])
+ >>> p.savez(outFile.name)
+ >>> _ = outFile.seek(0)
+ >>> p1 = tfields.Points3D.createFromFile(outFile.name)
+ >>> p.equal(p1)
+ True
+
+ """
+ kwargs = {}
+ for attr in self.__slots__:
+ if attr == '_cache':
+ # do not save the cache
+ continue
+ elif not hasattr(self, attr):
+ # attribute in __slots__ not found.
+ warnings.warn("When saving instance of class {0} Attribute {1} not set."
+ "This Attribute is not saved.".format(self.__class__, attr), Warning)
+ else:
+ kwargs[attr] = getattr(self, attr)
+ # kwargs = dict(zip(self.__slots__, [getattr(self, attr) for attr in self.__slots__]))
+
+ # resolve paths. Try except for buffer objects
+ try:
+ filePath = osTools.resolve(filePath)
+ except AttributeError:
+ pass
+ except:
+ raise
+ np.savez(filePath, self, **kwargs)
+
+ @classmethod
+ def createFromObjFile(cls, filePath, *groupNames):
+ """
+ Factory method
+ Given a filePath to a obj/wavefront file, construct the object
+ """
+ ioCls = ioTools.ObjFile
+ with ioCls(filePath, 'r') as f:
+ f.process()
+ vertices = f.getVertices(*groupNames)
+
+ log = logger.new()
+ if issubclass(cls, tfields.Triangles3D):
+ raise NotImplementedError("reading Triangles3D from obj. Use faces"
+ "attribute and the Mesh3D implementation.")
+ elif issubclass(cls, Points3D):
+ pass
+ else:
+ raise NotImplementedError("{cls} not implemented in createFromObjFile"
+ .format(**locals()))
+ return cls.__new__(cls, vertices)
+
+ @classmethod
+ def createFromTxtFile(cls, filePath):
+ """
+ Factory method
+ Given a filePath to a txt file, construct the object
+ """
+ ioCls = ioTools.HaukeFile
+ with ioCls(filePath, 'r') as f:
+ f.process()
+ return f.getPoints3D()
+
+ @classmethod
+ def createFromNpzFile(cls, filePath, **kwargs):
+ """
+ Factory method
+ Given a filePath to a npz file, construct the object
+ """
+ npFile = np.load(filePath)
+ keys = npFile.keys()
+ array = npFile['arr_0']
+ additionalKwargs = {key: npFile[key] for key in keys if key not in ['arr_0', '_cache']}
+ if additionalKwargs is not None:
+ kwargs.update(additionalKwargs)
+ return cls.__new__(cls, array, **kwargs)
+
+ @classmethod
+ def createFromInpFile(cls, filePath):
+ """
+ Factory method
+ Given a filePath to a inp file, construct the object
+ """
+ import ioTools.transcoding
+ transcoding = ioTools.transcoding.getTranscoding('inp')
+ content = transcoding.read(filePath)
+ part = content['parts'][0]
+ vertices = np.array([part['x'], part['y'], part['z']]).T / 1000
+ indices = np.array(part['nodeIndex']) - 1
+ if not list(indices) == range(len(indices)):
+ raise ValueError("node index skipped")
+ return cls(vertices)
+
+ @classmethod
+ def createFromFile(cls, filePath, *args, **kwargs):
+ """
+ load a file as a tensors object.
+ Args:
+ filePath (str or buffer)
+ *args:
+ forwarded to extension specific method
+ **kwargs:
+ extension (str): only needed if filePath is buffer
+ ... remaining:forwarded to extension specific method
+ """
+ extension = kwargs.pop('extension', 'npz')
+ if not osTools.isBuffer(filePath):
+ filePath = osTools.resolve(filePath)
+ extension = ioTools.getExtension(filePath)
+
+ # create from buffer, if file is tar file
+ if osTools.isTarPath(filePath):
+ with ioTools.TarHolder() as th:
+ buf = th.read(filePath)
+ inst = cls.createFromFile(buf, extension=extension, *args, **kwargs)
+ return inst
+
+ try:
+ fun = getattr(cls, 'createFrom{e}File'.format(e=extension.capitalize()))
+ except:
+ raise NotImplementedError("Do not understand the format of file {0}".format(filePath))
+ return fun(filePath, *args, **kwargs)
+
+ def indices(self, point):
+ """
+ Returns:
+ list of int: indices of point occuring
+ """
+ indices = []
+ for i, p in enumerate(self):
+ if all(p == point):
+ indices.append(i)
+ return indices
+
+ def index(self, point):
+ """
+ Args:
+ point
+ Returns:
+ int: index of point occuring
+ """
+ indices = self.indices(point)
+ if len(indices) == 0:
+ return None
+ if len(indices) == 1:
+ return indices[0]
+ raise ValueError("Multiple occurences of value %s" % point)
+
+ def mirrorCoordinate(self, coordinate, condition=None):
+ """
+ Examples:
+ >>> p = tfields.Points3D([[1., 2., 3.], [4., 5., 6.], [1, 2, -6]])
+ >>> p.mirrorCoordinate(1)
+ >>> p
+ Points3D([[ 1., -2., 3.],
+ [ 4., -5., 6.],
+ [ 1., -2., -6.]])
+
+ multiple coordinates can be mirrored. Eg. a point reflection would be
+ >>> p = tfields.Points3D([[1., 2., 3.], [4., 5., 6.], [1, 2, -6]])
+ >>> p.mirrorCoordinate([0,2])
+ >>> p
+ Points3D([[-1., 2., -3.],
+ [-4., 5., -6.],
+ [-1., 2., 6.]])
+
+ You can give a condition as mask or as str.
+ The mirroring will only be applied to the points meeting the condition.
+ >>> x, y, z = sympy.symbols('x y z')
+ >>> p.mirrorCoordinate([0,2], y > 3)
+ >>> p
+ Points3D([[-1., 2., -3.],
+ [ 4., 5., 6.],
+ [-1., 2., 6.]])
+
+ """
+ if condition is None:
+ condition = np.array([True for i in range(len(self))])
+ elif isinstance(condition, sympy.Basic):
+ condition = self.getMask(condition)
+ if isinstance(coordinate, list) or isinstance(coordinate, tuple):
+ for c in coordinate:
+ self.mirrorCoordinate(c, condition)
+ elif isinstance(coordinate, int):
+ self[:, coordinate][condition] *= -1
+ else:
+ raise TypeError()
+
+ def getMask(self, cutExpression=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.
+ Examples:
+ >>> 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]])
+ >>> p.getMask(x > 0)
+ array([ True, True, True, False, True, False], dtype=bool)
+
+ And combination
+ >>> p.getMask((x > 0) & (y < 3))
+ array([ True, False, True, False, True, False], dtype=bool)
+
+ Or combination
+ >>> p.getMask((x > 0) | (y > 3))
+ array([ True, True, True, False, True, False], dtype=bool)
+
+ """
+ coords = sympy.symbols('x y z')
+ with self.tmp_transform(coordSys):
+ mask = tfields.getMask(self, cutExpression, coords=coords)
+ return mask
+
+ def cut(self, cutExpression, coordSys=None):
+ """
+ Default cut method for Points3D. Works on a copy.
+ Args:
+ cutExpression (sympy logical expression): logical expression which will be evaluated.
+ use symbols x, y and z
+ coordSys (str): coordSys to evaluate the expression in.
+ Examples:
+ >>> 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]])
+ >>> p.cut(x > 0)
+ Points3D([[ 1., 2., 3.],
+ [ 4., 5., 6.],
+ [ 1., 2., -6.],
+ [ 1., 0., -1.]])
+
+ combinations of cuts
+ >>> p.cut((x > 0) & (z < 0))
+ Points3D([[ 1., 2., -6.],
+ [ 1., 0., -1.]])
+
+ Returns:
+ copy of self with cut applied
+
+ """
+ if len(self) == 0:
+ return self.copy()
+ mask = self.getMask(cutExpression, coordSys=coordSys)
+ mask.astype(bool)
+ inst = self[mask].copy()
+ return inst
+
+ def mapToCut(self, cutOrMeshMap, coordSys=None, atIntersection="remove"):
+ """
+ Args:
+ cutOrMeshMap (cutExpression / Mesh3D with faceScalars)
+ """
+ # in current implementation Points3D use always cut method, never MeshMap
+ return self.cut(cutOrMeshMap, coordSys=coordSys), None
+
+ def to3DArrays(self, arrayShape):
+ """
+ Args:
+ arrayShape (tuple of lenght 3): how coordinates will be devided.
+ If you initialize Points3D with a grid of dimensions a,b,c and
+ want to retrieve the grid-points use (a,b,c)
+ """
+ if not isinstance(arrayShape, tuple):
+ raise TypeError("ArrayShape is not of type tuple but %s" % type(arrayShape))
+ if not len(arrayShape) == 3:
+ raise TypeError("ArrayShape is not of length 3 but %i" % len(arrayShape))
+ return self.T.reshape((3,) + arrayShape)
+
+ def closestPoints(self, other, **kwargs):
+ """
+ Args:
+ other (Points3D): closest points to what? -> other
+ **kwargs: forwarded to scipy.spatial.cKDTree.query
+ Returns:
+ array shape(len(self)): Indices of other points that are closest to own points
+ Examples:
+ >>> m = tfields.Points3D([[1,0,0], [0,1,0], [1,1,0], [0,0,1],
+ ... [1,0,1]])
+ >>> p = tfields.Points3D([[1.1,1,0], [0,0.1,1], [1,0,1.1]])
+ >>> p.closestPoints(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)
+ array = res[1]
+
+ return array
+
+ def distances(self, other, metric='euclidean', **kwargs):
+ """
+ Examples:
+ >>> p = tfields.Points3D.createMeshGrid(np.linspace(0, 2, 3),
+ ... np.linspace(0, 2, 3),
+ ... np.linspace(0,0,1))
+ >>> p[4,2] = 1
+ >>> p.distances(p)[0,0]
+ 0.0
+ >>> p.distances(p)[5,1]
+ 1.4142135623730951
+ >>> p.distances([[0,1,2]])[-1]
+ array([ 3.])
+
+ """
+ return sp.spatial.distance.cdist(self, other, metric=metric, **kwargs)
+
+ def minDistances(self, other=None, metric='euclidean', **kwargs):
+ """
+ Examples:
+ >>> p = tfields.Points3D.createMeshGrid(np.linspace(0, 2, 3),
+ ... np.linspace(0, 2, 3),
+ ... np.linspace(0,0,1))
+ >>> p[4,2] = 1
+ >>> dMin = p.minDistances()
+ >>> dMin
+ array([ 1. , 1. , 1. , 1. , 1.41421356,
+ 1. , 1. , 1. , 1. ])
+
+ >>> dMin2 = p.minDistances(memorySaving=True)
+ >>> bool((dMin2 == dMin).all())
+ True
+
+ """
+ memorySaving = kwargs.pop('memorySaving', False)
+
+ if other is not None:
+ raise NotImplementedError("Should be easy but make shure not to remove 0s")
+ else:
+ other = self
+
+ try:
+ if memorySaving:
+ raise MemoryError()
+ d = self.distances(other, metric=metric, **kwargs)
+ return d[d > 0].reshape(d.shape[0], - 1).min(axis=1)
+ except MemoryError:
+ minDists = np.empty(self.shape[0])
+ for i, point in enumerate(other):
+ d = self.distances([point], metric=metric, **kwargs)
+ minDists[i] = d[d > 0].reshape(-1).min()
+ return minDists
+
+ def epsilonNeighbourhood(self, epsilon):
+ """
+ Returns:
+ indices for those sets of points that lie within epsilon around the other
+ Examples:
+ Create mesh grid with one extra point that will have 8 neighbours
+ within epsilon
+ >>> p = tfields.Points3D.createMeshGrid(np.linspace(0, 1, 2),
+ ... np.linspace(0, 1, 2),
+ ... np.linspace(0, 1, 2))
+ >>> p = tfields.Points3D([p, [0.5, 0.5, 0.5]])
+ >>> [len(en) for en in p.epsilonNeighbourhood(0.9)]
+ [2, 2, 2, 2, 2, 2, 2, 2, 9]
+
+ """
+ indices = np.arange(self.shape[0])
+ dists = self.distances(self)
+ distsInEpsilon = dists <= epsilon
+ return [indices[die] for die in distsInEpsilon]
+
+ def plot(self, **kwargs):
+ """
+ Frowarding to mpt.plotArray
+ """
+ artist = mpt.plotArray(self, **kwargs)
+ return artist
+
+
+if __name__ == '__main__':
+ import doctest
+ # doctest.testmod()
+ doctest.run_docstring_examples(Points3D, globals())
--
GitLab