From 51af578d997c2cb1e167fef5438a97bbacb7ccbc Mon Sep 17 00:00:00 2001 From: "Boeckenhoff, Daniel (dboe)" <daniel.boeckenhoff@ipp.mpg.de> Date: Tue, 14 Aug 2018 13:17:52 +0200 Subject: [PATCH] set_style generic --- test/test_core.py | 28 +++-- tfields/core.py | 26 ++++- tfields/mesh3D.py | 211 ++++++++++++++++++++++++++++++++--- tfields/plotting/__init__.py | 15 ++- tfields/plotting/mpl.py | 47 +++++--- tfields/triangles3D.py | 53 ++++++--- 6 files changed, 313 insertions(+), 67 deletions(-) diff --git a/test/test_core.py b/test/test_core.py index 2684a08..4e5b824 100644 --- a/test/test_core.py +++ b/test/test_core.py @@ -70,6 +70,19 @@ class Base_Test(object): def tearDown(self): del self._inst + + +class Tensor_Fields_Test(object): + def test_fields(self): + # field is of type list + self.assertTrue(isinstance(self._inst.fields, list)) + self.assertTrue(len(self._inst.fields) == len(self._fields)) + + for field, target_field in zip(self._inst.fields, self._fields): + self.assertTrue(np.array_equal(field, target_field)) + # fields are copied not reffered by a pointer + self.assertFalse(field is target_field) + """ EMPTY TESTS """ @@ -80,20 +93,10 @@ class Tensors_Empty_Test(Base_Test, unittest.TestCase): self._inst = tfields.Tensors([], dim=3) -class TensorFields_Empty_Test(Tensors_Empty_Test): +class TensorFields_Empty_Test(Tensors_Empty_Test, Tensor_Fields_Test): def setUp(self): self._fields = [] self._inst = tfields.TensorFields([], dim=3) - - def test_fields(self): - # field is of type list - self.assertTrue(isinstance(self._inst.fields, list)) - self.assertTrue(len(self._inst.fields) == len(self._fields)) - - for field, target_field in zip(self._inst.fields, self._fields): - self.assertTrue(np.array_equal(field, target_field)) - # fields are copied not reffered by a pointer - self.assertFalse(field is target_field) class TensorFields_Copy_Test(TensorFields_Empty_Test): @@ -104,6 +107,9 @@ class TensorFields_Copy_Test(TensorFields_Empty_Test): tensors = tfields.Tensors.grid(*base) self._inst = tfields.TensorFields(tensors, *self._fields) + self.assertTrue(self._fields[0].coord_sys, 'cylinder') + self.assertTrue(self._fields[1].coord_sys, 'cartesian') + class TensorMaps_Empty_Test(TensorFields_Empty_Test): def setUp(self): diff --git a/tfields/core.py b/tfields/core.py index 14d7daa..4ef3863 100644 --- a/tfields/core.py +++ b/tfields/core.py @@ -650,6 +650,9 @@ class Tensors(AbstractNdarray): ... [1, 4, 6], ... [1, 4, 7]]) True + + Lists or arrays are accepted also. + Furthermore, the iteration order can be changed >>> lins = tfields.Tensors.grid(np.linspace(3, 4, 2), np.linspace(0, 1, 2), ... np.linspace(6, 7, 2), iter_order=[1, 0, 2]) >>> lins.equal([[3, 0, 6], @@ -675,8 +678,22 @@ class Tensors(AbstractNdarray): ... [1, 4, 7]]) True - """ - inst = cls.__new__(cls, tfields.lib.grid.igrid(*base_vectors, **kwargs)) + When given the coord_sys argument, the grid is performed in the + given coorinate system: + >>> lins3 = tfields.Tensors.grid(np.linspace(4, 9, 2), + ... np.linspace(np.pi/2, np.pi/2, 1), + ... np.linspace(4, 4, 1), + ... iter_order=[2, 0, 1], + ... coord_sys=tfields.bases.CYLINDER) + >>> assert lins3.coord_sys == 'cylinder' + >>> lins3.transform('cartesian') + >>> assert np.array_equal(lins3[:, 1], [4, 9]) + + """ + cls_kwargs = {attr: kwargs.pop(attr) for attr in list(kwargs.keys()) if attr in cls.__slots__} + inst = cls.__new__(cls, + tfields.lib.grid.igrid(*base_vectors, **kwargs), + **cls_kwargs) return inst @property @@ -1144,9 +1161,10 @@ class Tensors(AbstractNdarray): """ indices = np.arange(self.shape[0]) - dists = self.distances(self) + dists = self.distances(self) # this takes long distsInEpsilon = dists <= epsilon - return [indices[die] for die in distsInEpsilon] + indices = [indices[die] for die in distsInEpsilon] # this takes long + return indices def _weights(self, weights, rigid=True): """ diff --git a/tfields/mesh3D.py b/tfields/mesh3D.py index 0237f03..e991267 100644 --- a/tfields/mesh3D.py +++ b/tfields/mesh3D.py @@ -9,7 +9,11 @@ part of tfields library import numpy as np import sympy import tfields + +# obj imports from tfields.lib.decorators import cached_property +import logging +import os def _dist_from_plane(point, plane): @@ -227,6 +231,148 @@ class Mesh3D(tfields.TensorMaps): raise ValueError("Face dimension should be 3") return obj + def _save_obj(self, path, **kwargs): + """ + Save obj as wavefront/.obj file + """ + obj = kwargs.pop('object', None) + group = kwargs.pop('group', None) + + cmap = kwargs.pop('cmap', 'viridis') + map_index = kwargs.pop('map_index', None) + + path = path.replace('.obj', '') + directory, name = os.path.split(path) + + if not (self.faceScalars.size == 0 or map_index is None): + import matplotlib.colors as colors + scalars = self.maps[0].fields[map_index] + min_scalar = scalars[~np.isnan(scalars)].min() + max_scalar = scalars[~np.isnan(scalars)].max() + vmin = kwargs.pop('vmin', min_scalar) + vmax = kwargs.pop('vmax', max_scalar) + if vmin == vmax: + if vmin == 0.: + vmax = 1. + else: + vmin = 0. + norm = colors.Normalize(vmin, vmax) + color_map = plt.get_cmap(cmap) + else: + # switch for not coloring the triangles and thus not producing the materials + norm = None + + if norm is not None: + mat_name = name + '_frame_{0}.mat'.format(map_index) + scalars[np.isnan(scalars)] = min_scalar - 1 + sorted_scalars = scalars[scalars.argsort()] + sorted_scalars[sorted_scalars == min_scalar - 1] = np.nan + sorted_faces = self.faces[scalars.argsort()] + scalar_set = np.unique(sorted_scalars) + scalar_set[scalar_set == min_scalar - 1] = np.nan + mat_path = os.path.join(directory, mat_name) + with open(mat_path, 'w') as mf: + for s in scalar_set: + if np.isnan(s): + mf.write("newmtl nan") + mf.write("Kd 0 0 0\n\n") + else: + mf.write("newmtl mtl_{0}\n".format(s)) + mf.write("Kd {c[0]} {c[1]} {c[2]}\n\n".format(c=color_map(norm(s)))) + else: + sorted_faces = self.faces + + # writing of the obj file + with open(path + '.obj', 'w') as f: + f.write("# File saved with tfields Mesh3D._save_obj method\n\n") + if norm is not None: + f.write("mtllib ./{0}\n\n".format(mat_name)) + 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)) + + last_scalar = None + for i, face in enumerate(sorted_faces + 1): + if norm is not None: + if not last_scalar == sorted_scalars[i]: + last_scalar = sorted_scalars[i] + f.write("usemtl mtl_{0}\n".format(last_scalar)) + f.write("f {f[0]} {f[1]} {f[2]}\n".format(f=face)) + + @classmethod + def _load_obj(cls, path, *group_names): + """ + Factory method + Given a path to a obj/wavefront file, construct the object + """ + import csv + log = logging.getLogger() + + with open(path, mode='r') as f: + reader = csv.reader(f, delimiter=' ') + groups = [] + group = None + vertex_no = 1 + for line in reader: + if not line: + continue + if line[0] == '#': + continue + if line[0] == 'g': + if group: + groups.append(group) + group = dict(name=line[1], vertices={}, faces=[]) + elif line[0] == 'v': + if not group: + log.warning("No group specified. I invent one myself.") + group = dict(name='Group', vertices={}, faces=[]) + vertex = list(map(float, line[1:4])) + group['vertices'][vertex_no] = vertex + vertex_no += 1 + elif line[0] == 'f': + face = [] + for v in line[1:]: + w = v.split('/') + face.append(int(w[0])) + group['faces'].append(face) + + vertices = [] + for g in groups[:]: + vertices.extend(g['vertices'].values()) + + if len(group_names) != 0: + groups = [g for g in groups if g['name'] in group_names] + + faces = [] + for g in groups: + faces.extend(g['faces']) + faces = np.add(np.array(faces), -1).tolist() + + """ + Building the class from retrieved vertices and faces + """ + if len(vertices) == 0: + return cls([]) + faceLenghts = [len(face) for face in faces] + for i in reversed(range(len(faceLenghts))): + length = faceLenghts[i] + if length == 3: + continue + if length == 4: + log.warning("Given a Rectangle. I will split it but " + "sometimes the order is different.") + faces.insert(i + 1, faces[i][2:] + faces[i][:1]) + faces[i] = faces[i][:3] + else: + raise NotImplementedError() + mesh = cls(vertices, faces=faces) + if group_names: + mesh = mesh.cleaned() + return mesh + @classmethod def plane(cls, *base_vectors, **kwargs): """ @@ -236,7 +382,7 @@ class Mesh3D(tfields.TensorMaps): be one-dimensional **kwargs: forwarded to __new__ """ - vertices = tfields.Tensors.grid(*base_vectors) + vertices = tfields.Tensors.grid(*base_vectors, **kwargs) base_vectors = tfields.grid.ensure_complex(*base_vectors) base_vectors = tfields.grid.to_base_vectors(*base_vectors) @@ -263,13 +409,41 @@ class Mesh3D(tfields.TensorMaps): idx_bot_right = len(base1) * (i0 + 1) + (i1 + 1) faces.append([idx_top_left, idx_top_right, idx_bot_left]) faces.append([idx_top_right, idx_bot_left, idx_bot_right]) - inst = cls.__new__(cls, vertices, faces=faces, **kwargs) + inst = cls.__new__(cls, vertices, faces=faces) return inst @classmethod def grid(cls, *base_vectors, **kwargs): """ Construct 'cuboid' along base_vectors + Examples: + Building symmetric geometries were never as easy: + + Approximated sphere with radius 1, translated in y by 2 units + >>> sphere = tfields.Mesh3D.grid((1, 1, 1), + ... (-np.pi, np.pi, 12), + ... (-np.pi / 2, np.pi / 2, 12), + ... coord_sys='spherical') + >>> sphere.transform('cartesian') + >>> sphere[:, 1] += 2 + + Oktaeder + >>> oktder = tfields.Mesh3D.grid((1, 1, 1), + ... (-np.pi, np.pi, 5), + ... (-np.pi / 2, np.pi / 2, 3), + ... coord_sys='spherical') + + Cube with edge length of 2 units + >>> cube = tfields.Mesh3D.grid((-1, 1, 2), + ... (-1, 1, 2), + ... (-5, -3, 2)) + + Cylinder + >>> cylinder = tfields.Mesh3D.grid((1, 1, 1), + ... (-np.pi, np.pi, 12), + ... (-5, 3, 12), + ... coord_sys='cylinder') + """ if not len(base_vectors) == 3: raise AttributeError("3 base_vectors vectors required") @@ -295,8 +469,7 @@ class Mesh3D(tfields.TensorMaps): basePart = base_vectors[:] basePart[coord] = np.array([base_vectors[coord][ind]], dtype=float) - - planes.append(cls.plane(*basePart)) + planes.append(cls.plane(*basePart, **kwargs)) inst = cls.merged(*planes, **kwargs) return inst @@ -352,16 +525,16 @@ class Mesh3D(tfields.TensorMaps): Check whether points lie within triangles with Barycentric Technique """ masks = self.triangles().in_triangles(points, delta, - assign_multiple=assign_multiple) + assign_multiple=assign_multiple) return masks - def removeFaces(self, faceDeleteMask): + def removeFaces(self, face_delete_mask): """ - Remove faces where faceDeleteMask is True + Remove faces where face_delete_mask is True """ - faceDeleteMask = np.array(faceDeleteMask, dtype=bool) - self.faces = self.faces[~faceDeleteMask] - self.faceScalars = self.faceScalars[~faceDeleteMask] + face_delete_mask = np.array(face_delete_mask, dtype=bool) + self.faces = self.faces[~face_delete_mask] + self.faceScalars = self.faceScalars[~face_delete_mask] def template(self, sub_mesh, delta=1e-9): """ @@ -384,13 +557,19 @@ class Mesh3D(tfields.TensorMaps): """ face_indices = np.arange(self.maps[0].shape[0]) cents = tfields.Tensors(sub_mesh.centroids()) - scalars = [] - mask = self.in_faces(cents, delta=delta) - scalars = [] - for face_mask in mask: - scalars.append(face_indices[face_mask][0]) + mask = self.in_faces(cents, delta) inst = sub_mesh.copy() - inst.maps[0].fields = [tfields.Tensors(scalars)] + if inst.maps: + scalars = [] + for face_mask in mask: + scalars.append(face_indices[face_mask][0]) + inst.maps[0].fields = [tfields.Tensors(scalars, dim=1)] + else: + inst.maps = [tfields.TensorFields([], + tfields.Tensors([], dim=1), + dim=3, + dtype=int) + ] return inst def _cut_sympy(self, expression, at_intersection="remove", _in_recursion=False): diff --git a/tfields/plotting/__init__.py b/tfields/plotting/__init__.py index 77d5d65..e6401ff 100644 --- a/tfields/plotting/__init__.py +++ b/tfields/plotting/__init__.py @@ -10,6 +10,7 @@ import matplotlib.pyplot as plt import matplotlib as mpl import numpy as np from .mpl import * +from six import string_types def set_default(dictionary, attr, value): @@ -121,7 +122,7 @@ class PlotOptions(object): vmax = self.get('vmax', vmaxDefault) return cmap, vmin, vmax - def formatColors(self, colors, fmt='rgba', length=None): + def format_colors(self, colors, fmt='rgba', length=None): """ format colors according to fmt argument Args: @@ -134,10 +135,11 @@ class PlotOptions(object): colors in fmt """ hasIter = True - if not hasattr(colors, '__iter__'): + if not hasattr(colors, '__iter__') or isinstance(colors, string_types): # colors is just one element hasIter = False colors = [colors] + if hasattr(colors[0], '__iter__') and fmt == 'norm': # rgba given but norm wanted cmap, vmin, vmax = self.getNormArgs(cmapDefault='NotSpecified', @@ -148,9 +150,9 @@ class PlotOptions(object): self.plotKwargs['vmax'] = vmax self.plotKwargs['cmap'] = cmap elif fmt == 'rgba': - if isinstance(colors[0], str) or isinstance(colors[0], unicode): + if isinstance(colors[0], string_types): # string color defined - colors = map(mpl.colors.to_rgba, colors) + colors = [mpl.colors.to_rgba(color) for color in colors] else: # norm given rgba wanted cmap, vmin, vmax = self.getNormArgs(cmapDefault='NotSpecified', @@ -218,3 +220,8 @@ class PlotOptions(object): if len(args) != 1: raise ValueError("Invalid number of args ({0})".format(len(args))) return self.retrieve(args[0], default, keep=keep) + + +if __name__ == '__main__': + import doctest + doctest.testmod() diff --git a/tfields/plotting/mpl.py b/tfields/plotting/mpl.py index 8588685..b3468a0 100644 --- a/tfields/plotting/mpl.py +++ b/tfields/plotting/mpl.py @@ -7,6 +7,7 @@ import numpy as np import warnings import os import matplotlib as mpl +from matplotlib import style import matplotlib.pyplot as plt from matplotlib.patches import Circle import mpl_toolkits.mplot3d as plt3D @@ -37,7 +38,7 @@ def gca(dim=None, **kwargs): return axis -def upgrade_style(style, source, dest="~/.config/matplotlib/"): +def upgrade_style(style, source, dest=None): """ Copy a style file at <origionalFilePath> to the <dest> which is the foreseen local matplotlib rc dir by default @@ -48,35 +49,42 @@ def upgrade_style(style, source, dest="~/.config/matplotlib/"): dest (str): local directory to copy the file to. Matpotlib has to search this directory for mplstyle files! """ - styleExtension = 'mplstyle' - path = tfields.lib.in_out.resolve(os.path.join(dest, style + '.' + styleExtension)) + if dest is None: + dest = mpl.get_configdir() + style_extension = 'mplstyle' + path = tfields.lib.in_out.resolve(os.path.join(dest, style + '.' + + style_extension)) source = tfields.lib.in_out.resolve(source) tfields.lib.in_out.cp(source, path) -def set_style(style='tfields', dest="~/.config/matplotlib/"): +def set_style(style='tfields', dest=None): """ Set the matplotlib style of name Important: Either you Args: style (str) - dest (str): local directory to use file from. if None, use default maplotlib styles + dest (str): local directory to use file from. if None, use default + maplotlib destination """ if dest is None: - path = style - else: - styleExtension = 'mplstyle' - path = tfields.lib.in_out.resolve(os.path.join(dest, style + '.' + styleExtension)) - try: + dest = mpl.get_configdir() + + style_extension = 'mplstyle' + path = tfields.lib.in_out.resolve(os.path.join(dest, style + '.' + + style_extension)) + if style in mpl.style.available: + plt.style.use(style) + elif os.path.exists(path): plt.style.use(path) - except IOError: + else: log = logging.getLogger() if style == 'tfields': log.warning("I will copy the default style to {dest}." .format(**locals())) source = os.path.join(os.path.dirname(__file__), - style + '.' + styleExtension) + style + '.' + style_extension) try: upgrade_style(style, source, dest) set_style(style) @@ -163,6 +171,7 @@ def plot_array(array, **kwargs): Artist or list of Artists (imitating the axis.scatter/plot behaviour). Better Artist not list of Artists """ + array = np.array(array) tfields.plotting.set_default(kwargs, 'methodName', 'scatter') po = tfields.plotting.PlotOptions(kwargs) @@ -194,6 +203,8 @@ def plot_mesh(vertices, faces, **kwargs): vmin vmax """ + vertices = np.array(vertices) + faces = np.array(faces) if faces.shape[0] == 0: warnings.warn("No faces to plot") return None @@ -241,9 +252,9 @@ def plot_mesh(vertices, faces, **kwargs): """ sort out color arguments """ - facecolors = po.formatColors(facecolors, - fmt='norm', - length=nFacesInitial) + facecolors = po.format_colors(facecolors, + fmt='norm', + length=nFacesInitial) if not full: facecolors = facecolors[dotProduct > 0] po.plotKwargs['facecolors'] = facecolors @@ -257,9 +268,9 @@ def plot_mesh(vertices, faces, **kwargs): color = po.retrieve_chain('color', 'c', 'facecolors', default='grey', keep=False) - color = po.formatColors(color, - fmt='rgba', - length=len(faces)) + color = po.format_colors(color, + fmt='rgba', + length=len(faces)) nanMask = np.isnan(color) if nanMask.any(): warnings.warn("nan found in colors. Removing the corresponding faces!") diff --git a/tfields/triangles3D.py b/tfields/triangles3D.py index d2a67af..c1cef3b 100644 --- a/tfields/triangles3D.py +++ b/tfields/triangles3D.py @@ -215,21 +215,21 @@ class Triangles3D(tfields.TensorFields): if not np.array_equal(transform, np.eye(3)): a = np.linalg.norm(np.linalg.solve(transform.T, - (self.bulk[aIndices, :] - - self.bulk[bIndices, :]).T), + (self[aIndices, :] - + self[bIndices, :]).T), axis=0) b = np.linalg.norm(np.linalg.solve(transform.T, - (self.bulk[aIndices, :] - - self.bulk[cIndices, :]).T), + (self[aIndices, :] - + self[cIndices, :]).T), axis=0) c = np.linalg.norm(np.linalg.solve(transform.T, - (self.bulk[bIndices, :] - - self.bulk[cIndices, :]).T), + (self[bIndices, :] - + self[cIndices, :]).T), axis=0) else: - a = np.linalg.norm(self.bulk[aIndices, :] - self.bulk[bIndices, :], axis=1) - b = np.linalg.norm(self.bulk[aIndices, :] - self.bulk[cIndices, :], axis=1) - c = np.linalg.norm(self.bulk[bIndices, :] - self.bulk[cIndices, :], axis=1) + a = np.linalg.norm(self[aIndices, :] - self[bIndices, :], axis=1) + b = np.linalg.norm(self[aIndices, :] - self[cIndices, :], axis=1) + c = np.linalg.norm(self[bIndices, :] - self[cIndices, :], axis=1) # sort by length for numerical stability lengths = np.concatenate((a.reshape(-1, 1), b.reshape(-1, 1), c.reshape(-1, 1)), axis=1) @@ -459,7 +459,9 @@ class Triangles3D(tfields.TensorFields): 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 + delta (float / None): + float: acceptance in +- norm vector direction + None: accept the face with the minimum distance to the point Returns: np.array: boolean mask, True where point in a triangle within delta Examples: @@ -487,6 +489,11 @@ class Triangles3D(tfields.TensorFields): ... m2.triangles()._in_triangles(np.array([0.2, 0.2, 0.1]), 0.2), ... np.array([ True, False], dtype=bool)) + if you set delta to None, the minimal distance point(s) are accepted + >>> assert np.array_equal( + ... m2.triangles()._in_triangles(np.array([0.2, 0.2, 0.1]), None), + ... np.array([ True, False], dtype=bool)) + If you define triangles that have colinear side vectors or in general lead to not invertable matrices the you will always get False >>> m3 = tfields.Mesh3D([[0,0,0], [2,0,0], [4,0,0], [0,1,0]], @@ -515,16 +522,33 @@ class Triangles3D(tfields.TensorFields): except: raise + # min_dist_method switch if delta is None + if delta is None: + delta = 1. + min_dist_method = True + else: + min_dist_method = False + u, v, w = self._baricentric(point, delta=delta).T if delta == 0.: w[np.isnan(w)] = 0. # division by 0 in baricentric makes w = 0 nan. + with warnings.catch_warnings(): warnings.filterwarnings('ignore', message="invalid value encountered in less_equal") - orthogonalAcceptance = (abs(w) <= 1) barycentric_bools = ((v <= 1.) & (v >= 0.)) & ((u <= 1.) & (u >= 0.)) & ((v + u <= 1.0)) - return np.array(barycentric_bools & orthogonalAcceptance) + if min_dist_method: + orthogonal_acceptance = np.full(barycentric_bools.shape, False) + closest_indices = np.argmin(abs(w)[barycentric_bools]) + # transform the indices to the whole array, not only the + # barycentric_bools selection + closest_indices = np.arange(len(barycentric_bools))[barycentric_bools][closest_indices] + orthogonal_acceptance[closest_indices] = True + else: + orthogonal_acceptance = (abs(w) <= 1) + + return np.array(barycentric_bools & orthogonal_acceptance) def in_triangles(self, tensors, delta=0., assign_multiple=False): """ @@ -562,6 +586,7 @@ class Triangles3D(tfields.TensorFields): masks = np.zeros((len(tensors), self.ntriangles()), dtype=bool) with tensors.tmp_transform(self.coord_sys): for i, point in enumerate(iter(tensors)): + # print(i) # print("Status: {i} points processed. Yet {nP} are found " # "to be in triangles (before handling assign_multiple)." # .format(nP=masks.sum(), **locals())) # SLOW! @@ -615,11 +640,11 @@ class Triangles3D(tfields.TensorFields): """ u, v, w = self._baricentric(point, 1.).T - orthogonalAcceptance = (w == 0) # point should lie in triangle + orthogonal_acceptance = (w == 0) # point should lie in triangle barycentric_bools = (((0. <= v) & (v <= 1.)) & (u == 0.)) | \ (((0. <= u) & (u <= 1.)) & (v == 0.)) | \ (v + u == 1.0) - return np.array(barycentric_bools & orthogonalAcceptance) + return np.array(barycentric_bools & orthogonal_acceptance) def _weights(self, weights, rigid=False): """ -- GitLab