From efadefd2d727da37d9560ce53a723a98b9723d2b Mon Sep 17 00:00:00 2001 From: Jan Caron <j.caron@fz-juelich.de> Date: Thu, 25 Jun 2020 20:02:09 +0200 Subject: [PATCH] Some Fixes and QoL Updates: vectors.py: More documentations and more fine control over vector and skyrmion chirality and orientation. Also fixed create_vector_singularity. colors: Added a namespace for colormaps (instead of a dictionary) to ease access to the maps. CMAP_CIRCULAR_DEFAULT is therefore deprecated now! plot2d: significantly simplified the symmetric color map handling with DivergingNorm (which was built for exactly that). --- src/empyre/fields/vectors.py | 96 ++++++++++++++++++++---------- src/empyre/io/field_plugins/vtk.py | 5 +- src/empyre/vis/colors.py | 35 +++++++---- src/empyre/vis/decorators.py | 6 +- src/empyre/vis/plot2d.py | 23 +++---- src/empyre/vis/plot3d.py | 11 +++- 6 files changed, 108 insertions(+), 68 deletions(-) diff --git a/src/empyre/fields/vectors.py b/src/empyre/fields/vectors.py index 9d4305a..a865a11 100644 --- a/src/empyre/fields/vectors.py +++ b/src/empyre/fields/vectors.py @@ -51,9 +51,7 @@ def create_vector_homog(dim, phi=0, theta=None, scale=1.0): return Field(data=data, scale=scale, vector=True) -def create_vector_vortex(dim, center=None, core_r=0, oop_r=None, axis=0, scale=1.0): - # TODO: oop and core_r documentation!!! General description! - # TODO: CHIRALITY +def create_vector_vortex(dim, center=None, phi_0=np.pi/2, oop_r=None, oop_sign=1, core_r=0, axis=0, scale=1.0): """Field subclass implementing a vortex vector field with 3 components in 2 or 3 dimensions. Attributes @@ -66,7 +64,31 @@ def create_vector_vortex(dim, center=None, core_r=0, oop_r=None, axis=0, scale=1 The vortex center should be between two pixels to avoid singularities. axis : int, optional The orientation of the vortex axis. The default is 0 and corresponds to the z-axis. Is ignored if dim is - only 2D. To invert chirality, multiply the resulting Field object by -1. # TODO: WRONG? + only 2D. + phi_0 : float, optional + Angular offset that allows all arrows to be rotated simultaneously. The default value is `pi/2` and corresponds + to an anti-clockwise rotation, while a value of `-pi/2` corresponds to a clockwise rotation. `0` and `pi` lead + to arrows that all point in/outward. + oop_r : float, optional + Radius of a potential out-of-plane ("oop") component in the vortex center. If `None`, the vortex is purely + in-plane, which is the default. Set this to a number (in pixels) that determines the the radius in which the + oop-component is tilted into the plane. A `0` leads to an immediate/sharp tilt (not visible without setting a + core radius, see `core_r`), while setting this to the radius of your vortex disc, will let it tilt smoothly + until reaching the edge. + oop_sign : {1, -1}, optional + Only used if `oop_r` is not None (i.e. there is an out-of-plane component). Can be `1` (default) or `-1` and + determines if the core (if it exists) is aligned parallel (`+1`) or antiparallel (`-1`) to the chosen symmetry + axis. + core_r : float, optional + Radius of a potential vortex core that's homogeneously oriented out-of-plane. Is not used when `oop_r` is not + set and defaults to `0`, which means the vortex core is infinitely sharp. + scale: tuple of float + Scaling along the dimensions of the underlying data. + + Returns + ------- + field : `~.Field` object + The resulting vector field. """ _log.debug('Calling create_vector_vortex') @@ -89,15 +111,19 @@ def create_vector_vortex(dim, center=None, core_r=0, oop_r=None, axis=0, scale=1 # Create vortex plane (2D): coords_uv = np.indices(dim_uv) + 0.5 # 0.5 to get to pixel/voxel center! coords_uv = coords_uv - np.asarray(center, dtype=float)[:, None, None] # Shift by center! - phi = np.arctan2(coords_uv[0], coords_uv[1]) - np.pi / 2 + phi = np.arctan2(coords_uv[0], coords_uv[1]) - phi_0 rr = np.hypot(coords_uv[0], coords_uv[1]) rr_clip = np.clip(rr - core_r, a_min=0, a_max=None) # negative inside core_r (clipped to 0), positive outside! + # Check if a core should be constructed (oop_r != 0): if oop_r is None: w_comp = np.zeros(dim_uv) else: - w_comp = 1 - 2/np.pi * np.arcsin(np.tanh(np.pi*rr_clip/oop_r)) # orthogonal: 1 inside, towards 0 outside! - v_comp = np.ones(dim_uv) * np.sin(phi) * np.sqrt(1 - w_comp) - u_comp = np.ones(dim_uv) * np.cos(phi) * np.sqrt(1 - w_comp) + assert oop_r >= 0, 'oop_r has to be a positive number!' + assert oop_sign in (1, -1), 'Sign of the out-of-plane component has to be either +1 or -1!' + w_comp = 1 - 2/np.pi * np.arcsin(np.tanh(np.pi*rr_clip/(oop_r+1E-30))) # orthogonal: 1 inside, to 0 outside! + w_comp *= oop_sign * w_comp + v_comp = np.ones(dim_uv) * np.sin(phi) * np.sqrt(1 - np.abs(w_comp)) + u_comp = np.ones(dim_uv) * np.cos(phi) * np.sqrt(1 - np.abs(w_comp)) if len(dim) == 3: # Expand to 3D: w_comp = np.expand_dims(w_comp, axis=axis) v_comp = np.expand_dims(v_comp, axis=axis) @@ -125,7 +151,7 @@ def create_vector_vortex(dim, center=None, core_r=0, oop_r=None, axis=0, scale=1 return Field(data=data, scale=scale, vector=True) -def create_vector_skyrmion(dim, center=None, phi_0=0, skyrm_d=None, wall_d=None, axis=0, scale=1.0): +def create_vector_skyrmion(dim, center=None, phi_0=0, core_sign=1, skyrm_d=None, wall_d=None, axis=0, scale=1.0): """Create a 3-dimensional magnetic Bloch or Neel type skyrmion distribution. Parameters @@ -139,20 +165,24 @@ def create_vector_skyrmion(dim, center=None, phi_0=0, skyrm_d=None, wall_d=None, phi_0 : float, optional Angular offset switching between Neel type (0 [default] or pi) or Bloch type (+/- pi/2) skyrmions. + core_sign : {1, -1}, optional + Can be `1` (default) or `-1` and determines if the skyrmion core is aligned parallel (`+1`) or antiparallel + (`-1`) to the chosen symmetry axis. skyrm_d : float, optional Diameter of the skyrmion. Defaults to half of the smaller dimension perpendicular to the skyrmion axis if not specified. wall_d : float, optional Diameter of the domain wall of the skyrmion. Defaults to `skyrm_d / 4` if not specified. - axis : {'z', '-z', 'y', '-y', 'x', '-x'}, optional # TODO: NUMBERS!! See vortex - The orientation of the skyrmion axis. The default is 'z'. Negative values invert skyrmion - core direction. + axis : int, optional + The orientation of the vortex axis. The default is 0 and corresponds to the z-axis. Is ignored if dim is + only 2D. + scale: tuple of float + Scaling along the dimensions of the underlying data. Returns ------- - amplitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3) - The magnetic distribution as a tuple of the 3 components in - `x`-, `y`- and `z`-direction on the 3-dimensional grid. + field : `~.Field` object + The resulting vector field. Notes ----- @@ -161,9 +191,8 @@ def create_vector_skyrmion(dim, center=None, phi_0=0, skyrm_d=None, wall_d=None, coordinates (e.g. coordinate 1 lies between the first and the second pixel). Skyrmion wall width is dependant on exchange stiffness A [J/m] and anisotropy K [J/m³] - The out-of-plane magnetization at the domain wall can be described as: - Mz = -Ms * tanh(x/w) # TODO: Instead ROMER paper - w = sqrt(A/K) + The out-of-plane magnetization at the domain wall is described in a paper by Romming et al + (see DOI: 10.1103/PhysRevLett.114.177203s) """ @@ -176,6 +205,7 @@ def create_vector_skyrmion(dim, center=None, phi_0=0, skyrm_d=None, wall_d=None, _log.debug('Calling create_vector_skyrmion') assert len(dim) in (2, 3), 'Skyrmion can only be build in 2 or 3 dimensions!' + assert core_sign in (1, -1), 'Sign of the out-of-plane component has to be either +1 or -1!' # Find indices of the skyrmion plane axes: idx_uv = [0, 1, 2] if len(dim) == 3: # 3D: @@ -201,7 +231,7 @@ def create_vector_skyrmion(dim, center=None, phi_0=0, skyrm_d=None, wall_d=None, rr = np.hypot(coords_uv[0], coords_uv[1]) phi = np.arctan2(coords_uv[0], coords_uv[1]) - phi_0 theta = _theta(rr) - w_comp = np.cos(theta) + w_comp = core_sign * np.cos(theta) v_comp = np.sin(theta) * np.sin(phi) u_comp = np.sin(theta) * np.cos(phi) # Expansion to 3D if necessary and component shuffling: @@ -243,16 +273,13 @@ def create_vector_singularity(dim, center=None, scale=1.0): The source center, given in 2D `(v, u)` or 3D `(z, y, x)`, where the perpendicular axis is discarded. Is set to the center of the field of view if not specified. The source center has to be between two pixels. - axis : {'z', '-z', 'y', '-y', 'x', '-x'}, optional # TODO: NUMBERS! - The orientation of the source axis. The default is 'z'. Negative values invert the source - to a sink. # TODO: wat...? - # TODO: scale! + scale: tuple of float + Scaling along the dimensions of the underlying data. Returns ------- - amplitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3) - The magnetic distribution as a tuple of the 3 components in - `x`-, `y`- and `z`-direction on the 3-dimensional grid. + field : `~.Field` object + The resulting vector field. Notes ----- @@ -260,7 +287,10 @@ def create_vector_singularity(dim, center=None, scale=1.0): reside at coordinates with _.5 at the end), i.e. integer values should be used as center coordinates (e.g. coordinate 1 lies between the first and the second pixel). - """ # TODO: What does negating do here? Senke / Quelle (YES IT DOES!)? swell and sink? ISSUE! INVITE PEOPLE!!! + Per default, all arrows point outwards, negate the resulting `Field` object with a minus after creation to let + them point inwards. + + """ _log.debug('Calling create_vector_singularity') # Find default values: if center is None: @@ -268,10 +298,12 @@ def create_vector_singularity(dim, center=None, scale=1.0): assert len(dim) == len(center), f"Length of dim ({len(dim)}) and center ({len(center)}) don't match!" # Setup coordinates, shape is (c, z, y, x), if 3D, or (c, y, x), if 2D (c: components): coords = np.indices(dim) + 0.5 # 0.5 to get to pixel/voxel center! - center = np.asarray(center, dtype=float) - bc_shape = (len(dim,),) + (1,)*len(dim) # Shape for broadcasting, (1,1,1,3) for 3D, (1,1,2) for 2D! - coords = coords - center.reshape(bc_shape) # Shift by center (append 1s for broadcasting)! + bc_shape = (len(dim,),) + (1,)*len(dim) # Shape for broadcasting, (3,1,1,1) for 3D, (2,1,1) for 2D! + coords = coords - np.reshape(center, bc_shape) # Shift by center (append 1s for broadcasting)! rr = np.sqrt(np.sum([coords[i]**2 for i in range(len(dim))], axis=0)) - data = coords / (rr + 1E-30) # Normalise amplitude (keep direction), rr (z,y,x) is broadcasted to data (c,z,y,x)! - data = data.T # (c,z,y,x) -> (x,y,z,c) + coords /= rr + 1E-30 # Normalise amplitude (keep direction), rr (z,y,x) is broadcasted to data (c,z,y,x)! + # coords has components listed in wrong order (z, y, x) in the first dimension, we need (x, y, z) in the last: + coords = np.flip(coords, axis=0) + # Finally, the data has the coordinate axis at the end, not at the front: + data = np.moveaxis(coords, 0, -1) # (c,z,y,x) -> (z,y,x,c) return Field(data=data, scale=scale, vector=True) diff --git a/src/empyre/io/field_plugins/vtk.py b/src/empyre/io/field_plugins/vtk.py index 1f91b7d..3be9a4a 100644 --- a/src/empyre/io/field_plugins/vtk.py +++ b/src/empyre/io/field_plugins/vtk.py @@ -118,7 +118,10 @@ def writer(filename, field, **kwargs): # Calculate colors: x_mag, y_mag, z_mag = field.comp magvec = np.asarray((x_mag.data.ravel(), y_mag.data.ravel(), z_mag.data.ravel())) - rgb = colors.CMAP_CIRCULAR_DEFAULT.rgb_from_vector(magvec) + cmap = kwargs.pop('cmap', None) + if cmap is None: + cmap = colors.cmaps.cyclic_cubehelix + rgb = cmap.rgb_from_vector(magvec) point_colors = tvtk.UnsignedIntArray() point_colors.number_of_components = 3 point_colors.name = 'colors' diff --git a/src/empyre/vis/colors.py b/src/empyre/vis/colors.py index 344a41f..6e71ecc 100644 --- a/src/empyre/vis/colors.py +++ b/src/empyre/vis/colors.py @@ -5,8 +5,7 @@ """This module provides a number of custom colormaps, which also have capabilities for 3D plotting. If this is the case, the :class:`~.Colormap3D` colormap class is a parent class. In `cmaps`, a -number of specialised colormaps is available for convenience. If the default for circular colormaps -(used for 3D plotting) should be changed, set it via `CMAP_CIRCULAR_DEFAULT`. +number of specialised colormaps is available for convenience. For general questions about colors see: http://www.poynton.com/PDFs/GammaFAQ.pdf http://www.poynton.com/PDFs/ColorFAQ.pdf @@ -28,7 +27,7 @@ from .tools import use_style __all__ = ['Colormap3D', 'ColormapCubehelix', 'ColormapPerception', 'ColormapHLS', 'ColormapClassic', - 'ColormapTransparent', 'cmaps', 'CMAP_CIRCULAR_DEFAULT', 'interpolate_color'] + 'ColormapTransparent', 'cmaps', 'interpolate_color'] _log = logging.getLogger(__name__) @@ -498,14 +497,24 @@ def interpolate_color(fraction, start, end): return r1, r2, r3 -cmaps = {'cubehelix_standard': ColormapCubehelix(), - 'cubehelix_reverse': ColormapCubehelix(reverse=True), - 'cubehelix_circular': ColormapCubehelix(start=1, rot=1, minLight=0.5, maxLight=0.5, sat=2), - 'perception_circular': ColormapPerception(), - 'hls_circular': ColormapHLS(), - 'classic_circular': ColormapClassic(), - 'transparent_black': ColormapTransparent(0, 0, 0, [0, 1.]), - 'transparent_white': ColormapTransparent(1, 1, 1, [0, 1.]), - 'transparent_confidence': ColormapTransparent(0.2, 0.3, 0.2, [0.75, 0.])} +class CMapNamespace(object): -CMAP_CIRCULAR_DEFAULT = cmaps['cubehelix_circular'] + def __init__(self): + self.cubehelix = ColormapCubehelix() + self.cubehelix_r = ColormapCubehelix(reverse=True) + self.cyclic_cubehelix = ColormapCubehelix(start=1, rot=1, minLight=0.5, maxLight=0.5, sat=2) + self.cyclic_perception = ColormapPerception() + self.cyclic_hls = ColormapHLS() + self.cyclic_classic = ColormapClassic() + self.transparent_black = ColormapTransparent(0, 0, 0, [0, 1.]) + self.transparent_white = ColormapTransparent(1, 1, 1, [0, 1.]) + self.transparent_confidence = ColormapTransparent(0.2, 0.3, 0.2, [0.75, 0.]) + + def __getitem__(self, key): + return self.__dict__[key] + + def add_cmap_dict(self, **kwargs): + self.__dict__.update(kwargs) + + +cmaps = CMapNamespace() diff --git a/src/empyre/vis/decorators.py b/src/empyre/vis/decorators.py index aa30034..f61990e 100644 --- a/src/empyre/vis/decorators.py +++ b/src/empyre/vis/decorators.py @@ -94,8 +94,8 @@ def colorwheel(axis=None, cmap=None, ax_size='20%', loc='upper right', **kwargs) axis : :class:`~matplotlib.axes.Axes`, optional Axis to which the colorwheel is added, by default None, which will pick the last used axis via `gca`. cmap : str or `matplotlib.colors.Colormap`, optional - The Colormap that should be used for the colorwheel, defaults to `None`, which chooses a map specified by - `.colors.CMAP_CIRCULAR_DEFAULT`. Needs to be a :class:`~.colors.Colormap3D`. + The Colormap that should be used for the colorwheel, defaults to `None`, which chooses the + `.colors.cmaps.cyclic_cubehelix` colormap. Needs to be a :class:`~.colors.Colormap3D` to work correctly. ax_size : str or float, optional String or float determining the size of the inset axis used, defaults to `20%`. loc : str or pair of floats, optional @@ -117,7 +117,7 @@ def colorwheel(axis=None, cmap=None, ax_size='20%', loc='upper right', **kwargs) ins_axes = inset_axes(axis, width=ax_size, height=ax_size, loc=loc) ins_axes.axis('off') if cmap is None: - cmap = colors.CMAP_CIRCULAR_DEFAULT + cmap = colors.cmaps.cyclic_cubehelix plt.sca(axis) # Set focus back to parent axis! return cmap.plot_colorwheel(axis=ins_axes, **kwargs) diff --git a/src/empyre/vis/plot2d.py b/src/empyre/vis/plot2d.py index 3b9abf5..19563bb 100644 --- a/src/empyre/vis/plot2d.py +++ b/src/empyre/vis/plot2d.py @@ -9,7 +9,7 @@ import logging import numpy as np import matplotlib.pyplot as plt -from matplotlib.colors import LinearSegmentedColormap +from matplotlib.colors import DivergingNorm from PIL import Image from . import colors @@ -75,19 +75,7 @@ def imshow(field, axis=None, cmap=None, **kwargs): elif isinstance(cmap, str): # make sure we have a Colormap object (and not a string): cmap = plt.get_cmap(cmap) if cmap.name.replace('_r', '') in DIVERGING_CMAPS: # 'replace' also matches reverted cmaps! - # Symmetric colormap only has zero at symmetry point if mappable has symmetric bounds (from - to + limit)! - vmin = kwargs.get('vmin', squeezed_field.data.min()-1E-30) - vmax = kwargs.get('vmax', squeezed_field.data.max()+1E-30) - vmin, vmax = np.min([vmin, 0]), np.max([0, vmax]) # Ensure zero is present! - kwargs.setdefault('vmin', vmin) - kwargs.setdefault('vmax', vmax) - limit = np.max(np.abs([vmin, vmax])) - # Calculate the subset of colors for the range vmin to vmax (often not full range of 2*limit): - start = 0.5 + vmin/(2*limit) # 0 for symmetric bounds, >0: unused colors at lower end! - end = 0.5 + vmax/(2*limit) # 1 for symmetric bounds, <1: unused colors at upper end! - cmap_colors = cmap(np.linspace(start, end, 256)) # New color indices with symmetry color at 0.5! - # Use calculated colors to create custom (asymmetric) colormap with symmetry color (often white) at zero: - cmap = LinearSegmentedColormap.from_list('custom', cmap_colors) + kwargs.setdefault('norm', DivergingNorm(0)) # Diverging colormap should have zero at the symmetry point! # Set extent in data coordinates (left, right, bottom, top) to kwargs (if not set explicitely): dim_v, dim_u, s_v, s_u = *squeezed_field.dim, *squeezed_field.scale kwargs.setdefault('extent', (0, dim_u * s_u, 0, dim_v * s_v)) @@ -199,7 +187,10 @@ def colorvec(field, axis=None, **kwargs): y_comp = comp[1] z_comp = comp[2] if (squeezed_field.ncomp == 3) else np.zeros(squeezed_field.dim) # Calculate image with color encoded directions: - rgb = colors.CMAP_CIRCULAR_DEFAULT.rgb_from_vector(np.stack((x_comp, y_comp, z_comp), axis=0)) + cmap = kwargs.pop('cmap', None) + if cmap is None: + cmap = colors.cmaps.cyclic_cubehelix + rgb = cmap.rgb_from_vector(np.stack((x_comp, y_comp, z_comp), axis=0)) # Set extent in data coordinates (left, right, bottom, top) to kwargs (if not set explicitely): dim_v, dim_u, s_v, s_u = *squeezed_field.dim, *squeezed_field.scale kwargs.setdefault('extent', (0, dim_u * s_u, 0, dim_v * s_v)) @@ -345,7 +336,7 @@ def quiver(field, axis=None, color_angles=False, cmap=None, n_bin='auto', bin_wi if color_angles: # Color angles according to calculated RGB values (only with circular colormaps): _log.debug('Encoding angles') if cmap is None: - cmap = colors.CMAP_CIRCULAR_DEFAULT + cmap = colors.cmaps.cyclic_cubehelix rgb = cmap.rgb_from_vector(np.asarray((x_comp, y_comp, z_comp))) / 255 rgba = np.concatenate((rgb, amplitude[..., None]), axis=-1) kwargs.setdefault('color', rgba.reshape(-1, 4)) diff --git a/src/empyre/vis/plot3d.py b/src/empyre/vis/plot3d.py index 54cda15..4afc9d0 100644 --- a/src/empyre/vis/plot3d.py +++ b/src/empyre/vis/plot3d.py @@ -109,7 +109,7 @@ def mask3d(field, title='Mask', threshold=0, grid=True, labels=True, return cont -def quiver3d(field, title='Vector Field', limit=None, cmap='jet', mode='2darrow', +def quiver3d(field, title='Vector Field', limit=None, cmap=None, mode='2darrow', coloring='angle', ar_dens=1, opacity=1.0, grid=True, labels=True, orientation=True, size=(700, 750), new_fig=True, view='isometric', position=None, bgcolor=(0.5, 0.5, 0.5)): @@ -122,7 +122,8 @@ def quiver3d(field, title='Vector Field', limit=None, cmap='jet', mode='2darrow' limit : float, optional Plotlimit for the vector field arrow length used to scale the colormap. cmap : string, optional - String describing the colormap which is used for amplitude encoding (default is 'jet'). + String describing the colormap which is used for color encoding (uses `~.colors.cmaps.cyclic_cubehelix` if + left on the `None` default) or amplitude encoding (uses 'jet' if left on the `None` default). ar_dens: int, optional Number defining the arrow density which is plotted. A higher ar_dens number skips more arrows (a number of 2 plots every second arrow). Default is 1. @@ -165,13 +166,17 @@ def quiver3d(field, title='Vector Field', limit=None, cmap='jet', mode='2darrow' vecs = mlab.quiver3d(xxx, yyy, zzz, x_mag, y_mag, z_mag, mode=mode, opacity=opacity, scalars=np.arange(len(xxx)), line_width=2) vector = np.asarray((x_mag.ravel(), y_mag.ravel(), z_mag.ravel())) - rgb = colors.CMAP_CIRCULAR_DEFAULT.rgb_from_vector(vector) + if cmap is None: + cmap = colors.cmaps.cyclic_cubehelix + rgb = cmap.rgb_from_vector(vector) rgba = np.hstack((rgb, 255 * np.ones((len(xxx), 1), dtype=np.uint8))) vecs.glyph.color_mode = 'color_by_scalar' vecs.module_manager.scalar_lut_manager.lut.table = rgba mlab.draw() elif coloring == 'amplitude': # Encodes the amplitude of the arrows with the jet colormap: _log.debug('Encoding amplitude') + if cmap is None: + cmap = 'jet' vecs = mlab.quiver3d(xxx, yyy, zzz, x_mag, y_mag, z_mag, mode=mode, colormap=cmap, opacity=opacity, line_width=2) mlab.colorbar(label_fmt='%.2f') -- GitLab