setup.py

-#!python
-# coding=utf-8
-"""Setup for testing, building, distributing and installing the 'Pyramid'-package"""
-
-import os
-import re
-import subprocess
-import sys
-from distutils.command.build import build
-
-import numpy
-from Cython.Distutils import build_ext
-from setuptools import setup, find_packages
-
-
-DISTNAME = 'pyramid'
-DESCRIPTION = 'PYthon based Reconstruction Algorithm for MagnetIc Distributions'
-MAINTAINER = 'Jan Caron'
-MAINTAINER_EMAIL = 'j.caron@fz-juelich.de'
-URL = ''
-VERSION = '0.1.0-dev'
-PYTHON_VERSION = (2, 7)
-DEPENDENCIES = {'numpy': (1, 10), 'cython': (0, 23)}
-LONG_DESCRIPTION = 'long description (TODO!)'  # TODO: Long description!
-
-
-def get_package_version(package):
-    """Return the package version of the specified package.
-
-    Parameters
-    ----------
-    package: basestring
-        Name of the package whic should be checked.
-
-    Returns
-    -------
-    version: tuple (N=3)
-        Version number as a tuple.
-
-    """
-    version = []
-    for version_attr in ('version', 'VERSION', '__version__'):
-        if (hasattr(package, version_attr) and
-                isinstance(getattr(package, version_attr), str)):
-            version_info = getattr(package, version_attr, '')
-            for part in re.split('\D+', version_info):
-                try:
-                    version.append(int(part))
-                except ValueError:
-                    pass
-    return tuple(version)
-
-
-def check_requirements():
-    """Checks the requirements of the Pyramid package."""
-    if sys.version_info < PYTHON_VERSION:
-        raise SystemExit('You need Python version %d.%d or later.'
-                         % PYTHON_VERSION)
-    for package_name, min_version in DEPENDENCIES.items():
-        dep_error = False
-        try:
-            package = __import__(package_name)
-        except ImportError:
-            dep_error = True
-        else:
-            package_version = get_package_version(package)
-            if min_version > package_version:
-                dep_error = True
-        if dep_error:
-            raise ImportError('You need `%s` version %d.%d or later.'
-                              % ((package_name,) + min_version))
-
-
-def hg_version():
-    """Get the Mercurial reference identifier.
-
-    Returns
-    -------
-    hg_ref: basestring
-        The Mercurial reference identifier.
-
-    """
-    try:
-        hg_rev = subprocess.check_output(['hg', 'id', '--id']).strip()
-    except:
-        hg_rev = "???"
-    return hg_rev
-
-
-def write_version_py(filename='pyramid/version.py'):
-    """Write the version.py file.
-
-    Parameters
-    ----------
-    filename: basestring, optional
-        Write the version and hg_revision into the specified python file.
-        Defaults to 'pyramid/version.py'.
-
-    """
-    version_string = '# -*- coding: utf-8 -*-\n' + \
-                     '""""This file is generated automatically by the Pyramid `setup.py`"""\n' + \
-                     'version = "{}"\n'.format(VERSION) + \
-                     'hg_revision = "{}"\n'.format(hg_version())
-    with open(os.path.join(os.path.dirname(__file__), filename), 'w') as vfile:
-        vfile.write(version_string)
-
-
-def get_files(rootdir):
-    """Returns a list of .py-files inside rootdir.
-
-    Parameters
-    ----------
-    rootdir: basestring
-        Root directory in which to search for ``.py``-files.
-
-    Returns
-    -------
-    filepaths: list
-        List of filepaths which were found.
-
-    """
-    filepaths = []
-    for root, dirs, files in os.walk(rootdir):
-        for filename in files:
-            if filename.endswith('.py'):
-                filepaths.append(os.path.join(root, filename))
-    return filepaths
-
-
-print('\n-------------------------------------------------------------------------------')
-print('checking requirements')
-check_requirements()
-print('write version.py')
-write_version_py()
-setup(name=DISTNAME,
-      description=DESCRIPTION,
-      long_description=LONG_DESCRIPTION,
-      maintainer=MAINTAINER,
-      maintainer_email=MAINTAINER_EMAIL,
-      url=URL,
-      download_url=URL,
-      version=VERSION,
-      packages=find_packages(exclude=['tests']),
-      include_dirs=[numpy.get_include()],
-      requires=['numpy', 'scipy', 'matplotlib', 'Pillow',
-                'mayavi', 'pyfftw', 'hyperspy', 'nose'],
-      test_suite='nose.collector',
-      cmdclass={'build_ext': build_ext, 'build': build})
-print('-------------------------------------------------------------------------------\n')

src/empyre/__init__.py

+# -*- coding: utf-8 -*-
+# Copyright 2020 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""Subpackage containing functionality for visualisation of multidimensional fields."""
+
+
+from . import fields
+from . import io
+from . import models
+from . import reconstruct
+from . import vis
+from . import utils
+
+__version__ = "0.3.4"
+__git_revision__ = "undefined"
+
+import logging
+_log = logging.getLogger(__name__)
+_log.info(f'Imported EMPyRe V-{__version__}')
+del logging
+
+__all__ = ['fields', 'io', 'models', 'reconstruct', 'vis', 'utils']

src/empyre/fields/__init__.py

+# -*- coding: utf-8 -*-
+# Copyright 2020 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""Subpackage containing container classes for multidimensional fields and ways to create them."""
+
+
+from .field import *
+from .shapes import *
+from .vectors import *
+
+
+__all__ = []
+__all__.extend(field.__all__)
+__all__.extend(shapes.__all__)
+__all__.extend(vectors.__all__)
+
+
+del field
+del shapes
+del vectors

src/empyre/fields/field.py

+# -*- coding: utf-8 -*-
+# Copyright 2020 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""This module provides a container class for multidimensional scalar or vector fields."""
+
+
+import logging
+from numbers import Number
+from numpy.lib.mixins import NDArrayOperatorsMixin
+
+import numpy as np
+try:
+    #  numpy >= 2
+    from numpy.lib.array_utils import normalize_axis_tuple
+except ModuleNotFoundError:
+    #  numpy < 2
+    from numpy.core.numeric import normalize_axis_tuple
+from scipy import ndimage
+
+from ..utils import Quaternion
+
+
+__all__ = ['Field']
+
+
+class Field(NDArrayOperatorsMixin):
+    """Container class for storing multidimensional scalar or vector fields.
+
+    The `Field` class is a sophisticated wrapper around a multidimensional numpy array. The user can access the
+    underlying data via the `data` attribute, or by using the `numpy.asarray` function.
+
+    `Field` defines the `ufunc` interface for numpys vast library of universal functions. This includes all operator
+    definitions, i.e., you can add, subtract, multiply and devide `Field` objects element-wise. Note that for vector
+    fields, `Field` will handle all vector components separately and for ufuncs with two inputs, the other object will
+    be broadcast accordingly. It is therefore also possible to add (or multiply) vector and scalar fields, assuming
+    their dimensions `dim` match (their `shape` can't match, because it includes the number of components `ncomp` for
+    the vector field). For ufuncs that reduce the output shape (e.g. `numpy.sum`), if the `axis` parameter is set to
+    `None` (default), component axes for vector fields are kept (they can however explicitely reduced by including
+    them in the `axis` paremeter, e.g. by setting `axis = -1`).
+
+    `Field` objects can be indexed like numpy arrays and scalar fields (`vector=False`) behave as expected. A vector
+    field (`vector=True`) can be indexed by only the non-component indices (see `dim`) and the according components of
+    the specified point are returned as a tuple. Slicing with integers will return a `Field` object with reduced
+    dimensionality (or a scalar/tuple if all spatial dimensions are reduced) that also drops the associated `scales`
+    as needed. New axes (indexing with `None` or `numpy.newaxis`) is disabled because that would mean that an unknown
+    scale would need to be added (you can always create a new `Field` object with an appropriately prepared ndarray).
+
+    Some callables exposed by numpy's public API (e.g. `numpy.mean(...)`) will treat `Field` objects as if they
+    were identical to their underlying numpy array (by envoking the `__array__` function which just returns `data`).
+    The result often is a simple `numpy.ndarray` and not a `Field` object. Some functions (e.g. `clip`) are also defined
+    on `Field` as methods. They might behave differently and will return `Field` objects instead, so take care and
+    refer to their docstrings for further information!
+    # TODO: Too verbose? In documentation instead?
+
+    Attributes
+    ----------
+    data: np.ndarray
+        The underlying data array of the field.
+    scale: tuple of float
+        Scaling along the dimensions of the underlying data.
+    vector: bool
+        True if the field is a vector field, False if it is a scalar field.
+    shape: tuple of ints
+        Shape of the underlying data. Includes the number of components as the first dimension if the field is a
+        vector field.
+    dim: tuple of ints
+        Dimensions of the underlying data. Only includes spatial dimensions, NOT the number of vector components!
+        (Use `ncomp` to get the number of components, if any, `shape` for the full shape of the underlying array).
+    ncomp: None or tuple of ints
+        Number of components if the field is a vector field (can be 2 or 3), None if it is a scalar field. The component
+        axis is always the last axis (index `-1`) of the underlying data array!
+
+    Notes
+    -----
+    See https://docs.scipy.org/doc/numpy/reference/arrays.classes.html for information about ndarray subclassing!
+    See https://docs.scipy.org/doc/numpy-1.13.0/neps/ufunc-overrides.html for information about __array_ufunc__!
+    See https://numpy.org/neps/nep-0018-array-function-protocol.html for information about __array_function__!
+    """
+
+    _log = logging.getLogger(__name__ + '.Field')
+
+    @property
+    def data(self):
+        return self.__data
+
+    @property
+    def vector(self):
+        return self.__vector
+
+    @vector.setter
+    def vector(self, vector):
+        assert isinstance(vector, bool), 'vector has to be a boolean!'
+        self.__vector = vector
+
+    @property
+    def scale(self):
+        return self.__scale
+
+    @scale.setter
+    def scale(self, scale):
+        if isinstance(scale, Number):  # Scale is the same for each dimension!
+            self.__scale = (float(scale),) * len(self.dim)
+        elif isinstance(scale, tuple):
+            ndim = len(self.dim)
+            assert len(scale) == ndim, f'Each of the {ndim} dimensions needs a scale, but {scale} was given!'
+            self.__scale = scale
+        else:
+            raise AssertionError(f'Scaling has to be a number or a tuple of numbers, was {scale} instead!')
+
+    @property
+    def shape(self):
+        return self.data.shape
+
+    @property
+    def dim(self):
+        if self.vector:
+            return self.shape[:-1]
+        else:
+            return self.shape
+
+    @property
+    def ncomp(self):
+        return self.shape[-1] if self.vector else 0
+
+    @property
+    def comp(self):  # TODO: Function or property? Philosophical question?
+        # Return empty list for scalars (ncomp is 0) and a list of components as scalar fields for a vector field:
+        return [Field(self.data[..., i], self.scale, vector=False) for i in range(self.ncomp)]
+
+    @property
+    def amp(self):  # TODO: Function or property? Philosophical question?
+        if self.vector:
+            amp = np.sqrt(np.sum([comp.data**2 for comp in self.comp], axis=0))
+        else:
+            amp = np.abs(self.data)
+        return Field(amp, self.scale, vector=False)
+
+    @property
+    def mask(self):  # TODO: Function or property? Philosophical question?
+        return Field(np.where(np.asarray(self.amp) > 0, True, False), self.scale, vector=False)
+
+    def __init__(self, data, scale=1.0, vector=False):
+        self._log.debug('Calling __init__')
+        self.__data = data
+        self.vector = vector  # Set vector before scale, because scale needs to know if vector (via calling dim)!
+        self.scale = scale
+        if self.vector:
+            assert self.ncomp in (2, 3), 'Only 2- or 3-component vector fields are supported!'
+        self._log.debug('Created ' + str(self))
+
+    def __repr__(self):
+        self._log.debug('Calling __repr__')
+        data_string = str(self.data)  # String of JUST the data array, without metadata (compared to repr)!
+        return f'{self.__class__.__name__}(data={data_string}, scale={self.scale}, vector={self.vector})'
+
+    def __str__(self):
+        self._log.debug('Calling __repr__')
+        ncomp_string = f', ncomp={self.ncomp}' if self.vector else ''
+        return f'{self.__class__.__name__}(dim={self.dim}, scale={self.scale}, vector={self.vector}{ncomp_string})'
+
+    def __array__(self):
+        self._log.debug('Calling __array__')
+        return self.data
+
+    def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
+        self._log.debug('Calling __array_ufunc__')
+        # outputs = kwargs.pop('out', ())  # TODO: Necessary?
+        outputs = kwargs.pop('out', (None,)*ufunc.nout)  # Defaults to tuple of None (currently: nout=1 all the time)
+        outputs_arr = tuple([np.asarray(out) if isinstance(out, Field) else out for out in outputs])
+        if all(a is None for a in outputs_arr):
+            outputs_arr = None
+        # Cannot handle items that have __array_ufunc__ (other than our own).
+        for item in inputs + outputs:
+            if hasattr(item, '__array_ufunc__') and not isinstance(item, Field):  # Something else with __array_ufunc__:
+                if type(item).__array_ufunc__ is not np.ndarray.__array_ufunc__:  # Can't handle other overrides:
+                    return NotImplemented
+        # TODO: BIGGEST NOTE HERE: Delegate work to ndarray.__array_ufunc__!
+        # TODO: For this to work, we have to make sure that we dispatch input as ndarrays, NOT Fields!
+        # Convert all Field inputs to simple ndarrays to avoid infinite recursion!
+        # ndarray.__array_ufunc__ delegates to the called ufunc method ONLY if all inputs and outputs are ndarrays
+        # (or have no __array_ufunc__ method, e.g. a scalar), so that's what we need to ensure here.:
+        # inputs = tuple([inp.view(np.ndarray) if isinstance(inp, Field) else inp for inp in inputs])
+        # TODO: possibly need to sort out which input (if two) is a Field (vector or scalar?), order matters?
+        # TODO: DOCUMENT: most things are same as in ndarrays, but if only one input is a vector field, try broadcast!
+        # TODO: for security, newaxis is not allowed (most other indexing works though), because scale would be unknown!
+        # 1 input (has to be a Field, otherwise we wouldn't be here):
+        self._log.debug(f'__array_ufunc__ inputs: {len(inputs)}')
+        self._log.info(f'ufunc: {ufunc}, method: {method}')
+        self._log.info(f'inputs: {inputs}')
+        self._log.info(f'outputs: {outputs}')
+        self._log.info(f'kwargs: {kwargs}')
+        if len(inputs) == 1:
+            field = inputs[0]
+            scale_new = field.scale
+            vector_new = field.vector
+            # Preprocess axis keyword if it exists:
+            axis = kwargs.get('axis', False)  # Default must not be None, because None is a possible setting!
+            full_reduction = False
+            # for ufunc.outer and ufunc.at:
+            if axis is not False:
+                ax_full = tuple(range(len(field.dim)))  # All axes (minus a possible component axis for vector Fields)!
+                ax_full_wc = tuple(range(len(field.dim) + 1))  # All axes WITH component axis (does not need to exist)!
+                axis = ax_full if axis is None else axis  # This keeps possible components untouched if axis was None!
+                # axis=-1 reduces over the vector components, if they exist
+                # Takes care of pot. neg. indices, ensures tuple!
+                axis = normalize_axis_tuple(axis, len(field.shape))
+                kwargs['axis'] = axis  # Put normalized axis back into kwargs!
+                if tuple(sorted(axis)) in (ax_full, ax_full_wc):
+                    full_reduction = True  # Full reduction (or reduction to just components) takes place:
+                if ax_full_wc[-1] in axis:  # User explicitely wants component reduction (can only be true for vector):
+                    vector_new = False  # Force scalar field!
+                scale_new = tuple([s for i, s in enumerate(field.scale) if i not in axis])  # Drop axis from scale!
+            inputs_arr = np.asarray(field)  # Convert inputs that are Fields to ndarrays to avoid recursion!
+            data_new = getattr(ufunc, method)(inputs_arr, out=outputs_arr, **kwargs)
+            if full_reduction:  # Premature return because the result is no longer a Field:
+                return data_new
+        # More than 1 input (at least one has to be a Field, otherwise we wouldn't be here):
+        elif len(inputs) > 1:
+            is_field = [isinstance(inp, Field) for inp in inputs]
+            is_vector = [getattr(inp, 'vector', False) for inp in inputs]
+            # Determine scale:
+            if np.sum(is_field) > 1:  # More than one input is a Field objects:
+                scales = [inp.scale for i, inp in enumerate(inputs) if is_field[i]]  # Only takes scales of Field obj.!
+                scale_new = scales[0]
+                err_msg = f'Scales of all Field objects must match! Given scales: {scales}!'
+                if not all(scale == scale_new for scale in scales):
+                    raise ValueError(err_msg)
+            else:  # Only one input is a field, pick the scale of that one:
+                scale_new = inputs[np.argmax(is_field)].scale  # argmax returns the index of first True!
+            # Determine vector:
+            vector_new = True if np.any(is_vector) else False  # Output is vector field if any input is a vector field!
+            if np.sum(is_vector) > 1:  # More than one input is a vector Field objects:
+                ncomps = [inp.ncomp for i, inp in enumerate(inputs) if is_vector[i]]  # Only takes ncomp of v.-Fields!
+                err_msg = f'# of components of all Field objects must match! Given ncomps: {ncomps}!'
+                if not all(ncomp == ncomps[0] for ncomp in ncomps):
+                    raise ValueError(err_msg)
+            # Append new axis at the end of non vector objects to broadcast to components:
+            if np.any(is_vector):
+                inputs = list(inputs)
+                for i, inp in enumerate(inputs):
+                    if not is_vector[i] and not isinstance(inp, Number):  # Numbers work for broadcasting anyway:
+                        if len(np.asarray(inp).shape) == len(scale_new):  # No. of dimensions w/o comp., have to match!
+                            inputs[i] = np.asarray(inputs[i])[..., np.newaxis]  # Broadcasting, try to cast as ndarray!
+                inputs = tuple(inputs)
+            # Convert inputs that are Fields to ndarrays to avoid recursion and determine data_new:
+            inputs_arr = tuple([np.asarray(inp) if isinstance(inp, Field) else inp for inp in inputs])
+            data_new = getattr(ufunc, method)(*inputs_arr, out=outputs_arr, **kwargs)
+        # Return results:
+        result = Field(data_new, scale_new, vector_new)
+        return result
+
+    def __getitem__(self, index):
+        self._log.debug('Calling __getitem__')
+        # Pre-processing index:
+        index = index if isinstance(index, tuple) else (index,)  # Make sure item is a tuple!
+        if None in index:
+            raise Exception('Field does not support indexing with newaxis/None! Please cast as ndarray first!')
+        if len(index) < len(self.shape) and Ellipsis not in index:  # Ellipsis is IMPLICIT at the end if index < dim:
+            index += (Ellipsis,)
+        if Ellipsis in index:  # Expand Ellipsis (...) into slice(None) (:) to make iterating consistent:
+            index = list(index)
+            i = index.index(Ellipsis)
+            missing_dims = len(self.shape) - len(index)  # Number of non-specified dimensions
+            index[i:i+1] = [slice(None)] * (missing_dims + 1)  # +1 because at least Ellipsis is replaced!
+            index = tuple(index)
+        # Indexing with a scalar drops the dimension, indexing with slice always keeps the dimension:
+        index_scale = index[:-1] if self.vector else index  # Disregard last index for vector fields (has no scale)!
+        scale_new = ()
+        for i, item in enumerate(index_scale):
+            if isinstance(item, slice):  # Indexing with slice keeps the dimension and therefore the scale:
+                scale_new += (self.scale[i],)
+        # Get data with classic indexing from underlying data array:
+        data_new = self.data[index]
+        # Determine vector state of output:
+        vector_new = self.vector
+        if self.vector and isinstance(index[-1], Number):  # For vector fields if last index (component) is dropped:
+            vector_new = False
+        # Return:
+        if not scale_new:  # scale_new=(), i.e. full reduction of all dimensions (not regarding possible vector comp.):
+            if isinstance(data_new, Number):  # full reduction:
+                return data_new
+            else:  # Only important for vector fields:
+                return tuple(data_new)  # return single vector as tuple
+        else:  # Return Field object
+            return Field(data_new, scale_new, vector=vector_new)
+
+    @classmethod
+    def from_scalar_fields(cls, scalar_list):
+        """Create a vector `Field` object from a list of scalar `Fields`.
+
+        Parameters
+        ----------
+        scalar_list : list
+            List of scalar `Field` objects (`vector=False`).
+
+        Returns
+        -------
+        vector_field: `Field`
+            A vector `Field` object containing the input scalar fields as components.
+
+        """
+        cls._log.debug('Calling from_scalar_fields')
+        assert len(scalar_list) in (2, 3), 'Needs two or three scalar fields as components for vector field!'
+        assert all(isinstance(item, Field) for item in scalar_list), 'All items have to be Field objects!'
+        assert all(not item.vector for item in scalar_list), 'All items have to be scalar fields!'
+        assert all(item.scale == scalar_list[0].scale for item in scalar_list), 'Scales of fields must match!'
+        return cls(np.stack(scalar_list, axis=-1), scalar_list[0].scale, vector=True)
+
+    @classmethod
+    def from_signal(cls, signal, scale=None, vector=None, comp_pos=-1):
+        """Convert a :class:`~hyperspy.signals.Signal` object to a :class:`~.Field` object.
+
+        Parameters
+        ----------
+        signal: :class:`~hyperspy.signals.Signal`
+            The :class:`~hyperspy.signals.Signal` object which should be converted to Field.
+
+        Returns
+        -------
+        field: :class:`~.Field`
+            A :class:`~.Field` object containing the loaded data.
+        scale: tuple of float, optional
+            Scaling along the dimensions of the underlying data. If not provided, will be read from the axes_manager.
+        vector: boolean, optional
+            If set to True, forces the signal to be interpreted as a vector field. EMPyRe will check if the first axis
+            is named 'vector components' (EMPyRe saves vector fields like this). If this is the case, vector will be
+            automatically set to True and the signal will also be interpreted as a vector field.
+        comp_pos: int, optoinal
+            The index of the axis containing the vector components (if `vector=True`). EMPyRe needs this to be the last
+            axis (index `-1`, which is the default). In case another position is given, the vector component will be
+            moved to the last axis. Old Pyramid files will have this axis at index `0`, so use this for backwards
+            compatibilty.
+
+        Notes
+        -----
+        Signals recquire the hyperspy package!
+
+        """
+        cls._log.debug('Calling from_signal')
+        data = signal.data
+        if vector and comp_pos != -1:  # component axis should be last, but is currently first -> roll to the end:
+            data = np.moveaxis(data, source=comp_pos, destination=-1)
+        if vector is None:  # Automatic detection:
+            vector = True if signal.axes_manager[0].name == 'vector components' else False
+        if scale is None:  # If not provided, try to read from axes_manager, one less axis if vector!:
+            scale = tuple([signal.axes_manager[i].scale for i in range(len(data.shape) - vector)])
+        return cls(data, scale, vector)
+
+    def to_signal(self):
+        """Convert :class:`~.Field` data into a HyperSpy signal.
+
+        Returns
+        -------
+        signal: :class:`~hyperspy.signals.Signal`
+            Representation of the :class:`~.Field` object as a HyperSpy Signal.
+
+        Notes
+        -----
+        This method recquires the hyperspy package!
+
+        """
+        self._log.debug('Calling to_signal')
+        try:  # Try importing HyperSpy:
+            import hyperspy.api as hs
+        except ImportError:
+            self._log.error('This method recquires the hyperspy package!')
+            return
+        # Create signal:
+        signal = hs.signals.BaseSignal(self.data)  # All axes are signal axes!
+        # Set axes:
+        if self.vector:
+            signal.axes_manager[0].name = 'vector components'
+        for i in range(len(self.dim)):
+            ax = i+1 if self.vector else i  # take component axis into account if vector!
+            num = ['x', 'y', 'z'][i] if len(self.dim) <= 3 else i
+            signal.axes_manager[ax].name = f'axis {num}'
+            signal.axes_manager[ax].scale = self.scale[i]
+            signal.axes_manager[ax].units = 'nm'
+        signal.metadata.Signal.title = f"EMPyRe {'vector' if self.vector else 'scalar'} Field"
+        return signal
+
+    def copy(self):
+        """Returns a copy of the :class:`~.Field` object.
+
+        Returns
+        -------
+        field: :class:`~.Field`
+            A copy of the :class:`~.Field`.
+
+        """
+        self._log.debug('Calling copy')
+        return Field(self.data.copy(), self.scale, self.vector)
+
+    def rotate(self, angle, axis='z', **kwargs):
+        """Rotate the :class:`~.Field`, based on :meth:`~scipy.ndimage.rotate`.
+
+        Rotation direction is from the first towards the second axis. Works for 2D and 3D Fields.
+
+        Parameters
+        ----------
+        angle : float
+            The rotation angle in degrees.
+        axis: {'x', 'y', 'z'}, optional
+            The axis around which the vector field is rotated. Default is 'z'. Ignored for 2D Fields.
+
+        Returns
+        -------
+        field: :class:`~.Field`
+            The rotated :class:`~.Field`.
+
+        Notes
+        -----
+        All additional kwargs are passed through to :meth:`~scipy.ndimage.rotate`.
+        The `reshape` parameter, controlling if the output shape is adapted so that the input array is contained
+        completely in the output is False per default, contrary to :meth:`~scipy.ndimage.rotate`,
+        where it is True.
+
+        """
+        self._log.debug('Calling rotate')
+        assert len(self.dim) in (2, 3), 'rotate is currently  only defined for 2D and 3D Fields!'
+        kwargs.setdefault('reshape', False)  # Default here is no reshaping!
+        if len(self.dim) == 2:  # For 2D, there are only 2 axes:
+            axis = 'z'  # Overwrite potential argument if 2D!
+            axes = (0, 1)  # y and x
+        else:  # 3D case:
+            # order of axes is important for permutation, to get positive levi_civita
+            axes = {'x': (0, 1), 'y': (2, 0), 'z': (1, 2)}[axis]
+        sc_0, sc_1 = self.scale[axes[0]], self.scale[axes[1]]
+        assert sc_0 == sc_1, f'rotate needs the scales in the rotation plane to be equal (they are {sc_0} & {sc_1})!'
+        np_angle = angle
+        if axis in ('x', 'z'):
+            np_angle *= -1
+        if not self.vector:  # Scalar field:
+            data_new = ndimage.rotate(self.data, np_angle, axes=axes, **kwargs)
+        else:  # Vector field:
+            # Rotate coordinate system:
+            comps = [np.asarray(comp) for comp in self.comp]
+            if self.ncomp == 3:
+                data_new = np.stack([ndimage.rotate(c, np_angle, axes=axes, **kwargs) for c in comps], axis=-1)
+                # Up till now, only the coordinates are rotated, now we need to rotate the vectors inside the voxels:
+                rot_axis = {'x': 2, 'y': 1, 'z': 0}[axis]
+                i, j, k = axes[0], axes[1], rot_axis  # next line only works if i != j != k
+                levi_civita = int((j-i) * (k-i) * (k-j) / (np.abs(j-i) * np.abs(k-i) * np.abs(k-j)))
+                # Create Quaternion, note that they have (x, y, z) order instead of (z, y, x):
+                quat_axis = tuple([levi_civita if i == rot_axis else 0 for i in (2, 1, 0)])
+                quat = Quaternion.from_axisangle(quat_axis, np.deg2rad(angle))
+                # T needed b.c. ordering!
+                data_new = quat.matrix.dot(data_new.reshape((-1, 3)).T).T.reshape(data_new.shape)
+            elif self.ncomp == 2:
+                u_comp, v_comp = comps
+                u_rot = ndimage.rotate(u_comp, np_angle, axes=axes, **kwargs)
+                v_rot = ndimage.rotate(v_comp, np_angle, axes=axes, **kwargs)
+                # Up till now, only the coordinates are rotated, now we need to rotate the vectors inside the voxels:
+                ang_rad = np.deg2rad(angle)
+                u_mix = np.cos(ang_rad)*u_rot - np.sin(ang_rad)*v_rot
+                v_mix = np.sin(ang_rad)*u_rot + np.cos(ang_rad)*v_rot
+                data_new = np.stack((u_mix, v_mix), axis=-1)
+            else:
+                raise ValueError('rotate currently only works for 2- or 3-component vector fields!')
+        return Field(data_new, self.scale, self.vector)
+
+    def rot90(self, k=1, axis='z'):
+        """Rotate the :class:`~.Field` 90° around the specified axis (right hand rotation).
+
+        Parameters
+        ----------
+        k : integer
+            Number of times the array is rotated by 90 degrees.
+        axis: {'x', 'y', 'z'}, optional
+            The axis around which the vector field is rotated. Default is 'z'. Ignored for 2D Fields.
+
+        Returns
+        -------
+        field: :class:`~.Field`
+            The rotated :class:`~.Field`.
+
+        """
+        self._log.debug('Calling rot90')
+        assert axis in ('x', 'y', 'z'), "Wrong input! 'x', 'y', 'z' allowed!"
+        assert len(self.dim) in (2, 3), 'rotate is currently  only defined for 2D and 3D Fields!'
+        if len(self.dim) == 2:  # For 2D, there are only 2 axes:
+            axis = 'z'  # Overwrite potential argument if 2D!
+            axes = (0, 1)  # y and x
+        else:  # 3D case:
+            axes = {'x': (0, 1), 'y': (0, 2), 'z': (1, 2)}[axis]
+        sc_0, sc_1 = self.scale[axes[0]], self.scale[axes[1]]
+        assert sc_0 == sc_1, f'rot90 needs the scales in the rotation plane to be equal (they are {sc_0} & {sc_1})!'
+        # TODO: Later on, rotation could also flip the scale (not implemented here, yet)!
+        if axis != 'y':
+            k = -k  # rotation is inverted if around x or z due to y flip compared to numpy
+        if not self.vector:  # Scalar Field:
+            data_new = np.rot90(self.data, k=k, axes=axes).copy()
+        else:  # Vector Field:
+            if len(self.dim) == 3:  # 3D:
+                assert self.ncomp == 3, 'rot90 currently only works for vector fields with 3 components in 3D!'
+                comp_x, comp_y, comp_z = self.comp
+                # reference:
+                # https://en.wikipedia.org/wiki/Rotation_matrix#In_three_dimensions
+                if axis == 'x':  # RotMatrix for 90°: [[1, 0, 0], [0, 0, -1], [0, 1, 0]]
+                    comp_x_rot = np.rot90(comp_x, k=k, axes=axes)
+                    comp_y_rot = -np.rot90(comp_z, k=k, axes=axes)
+                    comp_z_rot = np.rot90(comp_y, k=k, axes=axes)
+                elif axis == 'y':  # RotMatrix for 90°: [[0, 0, 1], [0, 1, 0], [-1, 0, 0]]
+                    comp_x_rot = np.rot90(comp_z, k=k, axes=axes)
+                    comp_y_rot = np.rot90(comp_y, k=k, axes=axes)
+                    comp_z_rot = -np.rot90(comp_x, k=k, axes=axes)
+                elif axis == 'z':  # RotMatrix for 90°: [[0, -1, 0], [1, 0, 0], [0, 0, 1]]
+                    comp_x_rot = -np.rot90(comp_y, k=k, axes=axes)
+                    comp_y_rot = np.rot90(comp_x, k=k, axes=axes)
+                    comp_z_rot = np.rot90(comp_z, k=k, axes=axes)
+                data_new = np.stack((comp_x_rot, comp_y_rot, comp_z_rot), axis=-1)
+            if len(self.dim) == 2:  # 2D:
+                assert self.ncomp == 2, 'rot90 currently only works for vector fields with 2 components in 2D!'
+                comp_x, comp_y = self.comp
+                comp_x_rot = -np.rot90(comp_y, k=k, axes=axes)
+                comp_y_rot = np.rot90(comp_x, k=k, axes=axes)
+                data_new = np.stack((comp_x_rot, comp_y_rot), axis=-1)
+        # Return result:
+        return Field(data_new, self.scale, self.vector)
+
+    def get_vector(self, mask=None):
+        """Returns the field as a vector, specified by a mask.
+
+        Parameters
+        ----------
+        mask : :class:`~numpy.ndarray` (boolean)
+            Masks the pixels from which the entries should be taken. Must be a numpy array for indexing to work!
+
+        Returns
+        -------
+        vector : :class:`~numpy.ndarray` (N=1)
+            The vector containing the field of the specified pixels. If the field is a vector field, components are
+            first vectorised, then concatenated!
+
+        """
+        self._log.debug('Calling get_vector')
+        if mask is None:  # If not given, take full data:
+            mask = np.full(self.dim, True)
+        if self.vector:  # Vector field:
+            return np.ravel([comp.data[mask] for comp in self.comp])
+        else:  # Scalar field:
+            return np.ravel(self.data[mask])
+
+    def set_vector(self, vector, mask=None):
+        """Set the field of the masked pixels to the values specified by `vector`.
+
+        Parameters
+        ----------
+        mask : :class:`~numpy.ndarray` (boolean), optional
+            Masks the pixels from which the field should be taken.
+        vector : :class:`~numpy.ndarray` (N=1)
+            The vector containing the field of the specified pixels.
+
+        Returns
+        -------
+        None
+
+        """
+        self._log.debug('Calling set_vector')
+        if mask is None:  # If not given, set full data:
+            mask = np.full(self.dim, True)
+        if self.vector:  # Vector field:
+            assert np.size(vector) % self.ncomp == 0, 'Vector has to contain all components for every pixel!'
+            count = np.size(vector) // self.ncomp
+            for i in range(self.ncomp):
+                sl = slice(i*count, (i+1)*count)
+                self.data[..., i][mask] = vector[sl]
+        else:  # Scalar field:
+            self.data[mask] = vector
+
+    def squeeze(self):
+        """Squeeze the `Field` object to remove single-dimensional entries in the shape.
+
+        Returns
+        -------
+        field: `Field`
+            Squeezed `Field` object. Note that scales associated with squeezed dimensions are also dropped.
+
+        Notes
+        -----
+        Also works via `numpy.squeeze(field)`, because `numpy.squeeze` searches for a local implementation first and
+        then uses `_wrapit` to envoke this function here!
+
+        """
+        self._log.debug('Calling squeeze')
+        squeezed_indices = np.flatnonzero(np.asarray(self.dim) == 1)
+        squeezed_data = np.squeeze(self.data)
+        if squeezed_indices:
+            self._log.info(f'The following indices were squeezed: {squeezed_indices}')
+        squeezed_scale = tuple([self.scale[i] for i in range(len(self.dim)) if i not in squeezed_indices])
+        return Field(squeezed_data, squeezed_scale, self.vector)
+
+    def pad(self, pad_width, mode='constant', **kwargs):
+        """Pad the `Field` object and increase the size of the underlying array.
+
+        Parameters
+        ----------
+        pad_width : {sequence, array_like, int}
+            Number of values padded to the edges of each axis. Can be a single number for all, one number per axis or
+            a tuple `(before, after)` for each axis for finer control.
+        mode : str, optional
+            Padding mode (see `numpy.pad` for all options), by default 'constant', which pads with zeros.
+
+        Returns
+        -------
+        field: `Field`
+            The padded `Field` object.
+
+        """
+        if isinstance(pad_width, Number):  # Paddding is the same for each dimension (make sure it is a tuple)!
+            pad_width = (pad_width,) * len(self.dim)
+        pad_width = [(p, p) if isinstance(p, Number) else p for p in pad_width]
+        if self.vector:  # Append zeros to padding, so component axis stays as is:
+            pad_width = pad_width + [(0, 0)]
+        data_new = np.pad(self.data, pad_width, mode, **kwargs)
+        return Field(data_new, self.scale, self.vector)
+
+    def bin(self, n):
+        """Bins data of the `Field` to decrease the size of the underlying array by averaging over a number of values.
+
+        Parameters
+        ----------
+        n : {sequence, array_like, int}
+            Number of entries along each axis over which the average is taken. Can be a single integer for all axes or a
+            list, specifying the number of entries over which to average for each individual axis.
+
+        Returns
+        -------
+        field: `Field`
+            The binned `Field` object.
+
+        Notes
+        -----
+        Padding takes place before binning to ensure the dimensions are multiples of `n`. The padding mode is `edge`.
+
+        """
+        self._log.debug('Calling bin')
+        assert isinstance(n, (int, tuple)), 'n must be a positive integer or a tuple of positive integers!'
+        if isinstance(n, Number):  # Binning is the same for each dimension (make sure it is a tuple)!
+            n = (n,) * len(self.dim)
+        assert all([n_i >= 0 for n_i in n]), 'All entries of n must be positive integers!'
+        # Pad if necessary (use padded 'field' from here on), formula for clarity: (n - dim % n) % n
+        pad_width = [(0, (n[i] - self.dim[i] % n[i]) % n[i]) for i in range(len(self.dim))]
+        field = self.pad(pad_width, mode='edge')
+        # Create new shape used for binning (mean over every second axis will be taken):
+        bin_shape = list(np.ravel([(field.dim[i]//n[i], n[i]) for i in range(len(field.dim))]))
+        mean_axes = np.arange(1, 2*len(field.dim), 2)  # every 2nd axis!
+        if self.vector:  # Vector field:
+            bin_shape += [field.ncomp]  # Append component axis (they stay unchanged)
+        # Bin data and scale accordingly:
+        data_new = field.data.reshape(bin_shape).mean(axis=tuple(mean_axes))
+        scale_new = tuple([field.scale[i] * n[i] for i in range(len(field.dim))])
+        return Field(data_new, scale_new, self.vector)
+
+    def zoom(self, zoom, **kwargs):
+        """Wrapper for the `scipy.ndimage.zoom` function.
+
+        Parameters
+        ----------
+        zoom : float or sequence
+            Zoom factor along the axes. If a float, `zoom` is the same for each axis. If a sequence, `zoom` needs to
+            contain one value for each axis.
+
+        Returns
+        -------
+        field: `Field`
+            The zoomed in `Field` object.
+
+        Notes
+        -----
+        As in `scipy.ndimage.zoom`, a spline order can be specified, which defaults to 3.
+
+        """
+        self._log.debug('Calling zoom')
+        if not self.vector:  # Scalar field:
+            data_new = ndimage.zoom(self.data, zoom, **kwargs)
+        else:  # Vector field:
+            comps = [np.asarray(comp) for comp in self.comp]
+            data_new = np.stack([ndimage.zoom(comp, zoom, **kwargs) for comp in comps], axis=-1)
+        if isinstance(zoom, Number):  # Zoom is the same for each dimension!
+            zoom = (zoom,) * len(self.dim)
+        scale_new = tuple([self.scale[i]/z for i, z in enumerate(zoom)])
+        return Field(data_new, scale_new, self.vector)
+
+    def flip(self, axis=None, **kwargs):
+        """Returns a `Field` with entries flipped along specified axes.
+
+        Parameters
+        ----------
+        axis : tuple or None, optional
+            Axes whose entries will be flipped, by default None.
+        """
+        self._log.debug('Calling flip')
+        if self.vector and axis is None:
+            axis = tuple(range(len(self.dim)))
+        axis = normalize_axis_tuple(axis, len(self.shape))  # Make sure, axis is a tuple!
+        data_new = np.flip(self.data, axis, **kwargs).copy()  # Only flips space, not component direction!
+        if self.vector:
+            flip_vec = [
+                -1 if i in axis else 1
+                for i in reversed(range(self.ncomp))
+            ]
+            data_new *= flip_vec
+        return Field(data_new, self.scale, self.vector)
+
+    def gradient(self):
+        """Returns the gradient of an N-dimensional scalar `Field`. Wrapper around the according numpy function.
+
+        Returns
+        -------
+        gradients: ndarray or list of ndarray
+            A set of ndarrays (or a single ndarray for 1D input), corresponding to the derivatives of the field with
+            respect to each dimension.
+
+        Notes
+        -----
+        The field is implicitely squeezed before the gradient is calculated!
+
+        """
+        self._log.debug('Calling gradient')
+        assert not self.vector, 'Gradient can only be computed from scalar fields!'
+        squeezed_field = self.squeeze()  # Gradient along dimension of length 1 does not work!
+        gradients = np.gradient(np.asarray(squeezed_field), *squeezed_field.scale)
+        if len(squeezed_field.dim) == 1:  # Only one gradient!
+            return Field(gradients, squeezed_field.scale, vector=False)
+        else:  # Multidimensional gradient (flip component order with [::-1], so that e.g x/y/z instead of z/y/x):
+            return Field(np.stack(gradients[::-1], axis=-1), squeezed_field.scale, vector=True)
+
+    def curl(self):
+        """Returns the curl of an N-dimensional `Field`.
+
+        Returns
+        -------
+        field: `Field`
+            The curl/rotation.
+
+        Notes
+        -----
+        The calculation depends on the input:
+        3 dimensions, 3 components: Calculates the full 3D rotational vector field!
+        2 dimensions, 2 components: Calculates the out-of-plane component of the curl as a 2D scalar field!
+        2 dimensions, scalar: Calculates the planar rotation as a 2D vector field!
+
+        """
+        self._log.debug('Calling curl')
+        squeezed_field = self.squeeze()
+        if len(squeezed_field.dim) == 3:  # 3 dimensions:
+            if squeezed_field.ncomp == 3:  # 3 component vector field (standard case):
+                self._log.debug('input: 3 dimensions, 3 components!')
+                field_x, field_y, field_z = squeezed_field.comp
+                gradx_x, grady_x, gradz_x = field_x.gradient().comp
+                gradx_y, grady_y, gradz_y = field_y.gradient().comp
+                gradx_z, grady_z, gradz_z = field_z.gradient().comp
+                curl_x = grady_z - gradz_y
+                curl_y = gradz_x - gradx_z
+                curl_z = gradx_y - grady_x
+                return Field.from_scalar_fields([curl_x, curl_y, curl_z])
+            else:
+                raise AssertionError('Can only handle 3 component vector fields in 3D!')
+        elif len(squeezed_field.dim) == 2:  # 2 dimensions (tricky, usually not hardly defined):
+            if squeezed_field.ncomp == 2:  # 2 component vector field (return perpendicular component as scalar field):
+                self._log.debug('input: 2 dimensions, 2 components!')
+                field_x, field_y = squeezed_field.comp
+                gradx_x, grady_x = field_x.gradient().comp
+                gradx_y, grady_y = field_y.gradient().comp
+                return gradx_y - grady_x
+            elif not squeezed_field.vector:  # scalar field (return planar components as 2D vector field):
+                self._log.debug('input: 3 dimensions, scalar field!')
+                gradx, grady = squeezed_field.gradient().comp
+                curl_x = grady
+                curl_y = -gradx
+                return Field.from_scalar_fields([curl_x, curl_y])
+            else:
+                raise AssertionError('Can only handle 3 component vector or scalar fields in 2D!')
+        else:
+            raise AssertionError('Can only handle 3D or 2D cases (see documentation for specifics)!')
+
+    def clip(self, vmin=None, vmax=None, sigma=None, mask=None):
+        """Clip (limit) the values in an array. For vector fields, this will take the amplitude into account.
+
+        Parameters
+        ----------
+        vmin : float, optional
+            Mimimum value, by default None. Is ignored for vector fields. Will overwrite values found via `sigma` or
+            `mask`.
+        vmax : float, optional
+            Maximum value, by default None. Will overwrite values found via `sigma` or `mask`.
+        sigma : float, optional
+            Defines a range in standard deviations, that results in a boolean mask that excludes outliers, by default
+            None. E.g. `sigma=2` will mark points within the 5% highest amplitude as outliers. `vmin` and `vmax` will be
+            searched in the remaining region. Is logically combined with `mask` if both are set.
+        mask : ndarray, optional
+            Boolean array that directly describes where to search for `vmin` and `vmax`, by default None. Is logically
+            combined with the mask from `sigma` if both are set.
+
+        Returns
+        -------
+        field: `Field`
+            The clipped `Field`.
+
+        Notes
+        -----
+        In contrast to the corresponding numpy function, `vmin` and `vmax` can both be `None` due to the alternative
+        clipping strategies employed here.
+
+        """
+        self._log.debug('Calling clip')
+        # Get a scalar indicator array for where clipping needs to happen:
+        if self.vector:  # For vector fields, it is the amplitude of the data:
+            indicator = self.amp.data
+        else:  # For scalar fields this is the data itself:
+            indicator = self.data
+        if mask is None:  # If no mask is set yet, default to True everywhere:
+            mask = np.full(self.dim, True)
+        if sigma is not None:  # Mark outliers that are outside `sigma` standard deviations:
+            sigma_mask = (indicator - indicator.mean()) < (sigma * indicator.std())
+            mask = np.logical_and(mask, sigma_mask)
+        # Determine vmin and vmax if they are not set by the user:
+        indicator_masked = np.where(mask, indicator, np.nan)
+        if vmin is None:
+            vmin = np.nanmin(indicator_masked)
+        if vmax is None:
+            vmax = np.nanmax(indicator_masked)
+        if self.vector:  # Vector fields need to scale components according to masked amplitude:
+            # Only vmax is important for vectors! mask_vec broadcast to components!
+            mask_vec = (indicator <= vmax)[..., None]
+            data_new = np.where(mask_vec, self.data, vmax * self.data/indicator[..., None])  # Scale outliers to vmax!
+        else:  # For scalar fields, just delegate to the numpy function:
+            data_new = np.clip(self.data, vmin, vmax)
+        return Field(data_new, self.scale, self.vector)

src/empyre/fields/shapes.py

+# -*- coding: utf-8 -*-
+# Copyright 2020 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+
+
+import logging
+
+import numpy as np
+
+from .field import Field
+
+
+__all__ = ['create_shape_slab', 'create_shape_disc', 'create_shape_ellipse', 'create_shape_ellipsoid',
+           'create_shape_sphere', 'create_shape_filament', 'create_shape_voxel']
+
+_log = logging.getLogger(__name__)
+
+
+# TODO: TEST!!!
+
+
+def create_shape_slab(dim, center=None, width=None, scale=1.0):
+    """Creates a Field object with the shape of a slab as a scalar field in arbitrary dimensions.
+
+    Attributes
+    ----------
+    dim : tuple
+        The dimensions of the grid.
+    center : tuple, optional
+        The center of the slab in pixel coordinates.
+    width : tuple, optional
+        The width of the slab in pixel coordinates.
+    scale: tuple of float
+        Scaling along the dimensions of the underlying data.
+
+    """
+    _log.debug('Calling __init__')
+    if center is None:
+        center = tuple([d/2 for d in dim])
+    if width is None:
+        width = tuple([d/2 for d in dim])
+    assert len(dim) == len(center) == len(width), 'Parameters dim, center and width must have the same dimensions!'
+    data = np.zeros(dim)
+    bounds = ()
+    for i in range(len(dim)):
+        start = int(np.floor(center[i] - width[i]/2))
+        stop = int(np.ceil(center[i] + width[i]/2))
+        bounds += (slice(start, stop),)
+    data[bounds] = 1
+    return Field(data=data, scale=scale, vector=False)
+
+
+def create_shape_disc(dim, center=None, radius=None, height=None, axis=0, scale=1.0):
+    """Creates a Field object with the shape of a cylindrical disc in 2D or 3D.
+
+    Attributes
+    ----------
+    dim : tuple
+        The dimensions of the grid.
+    center : tuple, optional
+        The center of the disc in pixel coordinates.
+    radius : float, optional
+        The radius of the disc in pixel coordinates.
+    height : float, optional
+        The height of the disc in pixel coordinates. Unused if only 2D.
+    axis : int, optional
+        The orientation of the discs orthogonal axis. Only used in 3D case with z-axis as default.
+    scale: tuple of float
+        Scaling along the dimensions of the underlying data.
+
+    """
+    _log.debug('Calling __init__')
+    assert len(dim) in (2, 3), 'Disc can only be build in 2 or 3 dimensions!'
+    # Find indices of the disc plane axes:
+    idx_uv = [0, 1, 2]
+    if len(dim) == 3:  # 3D:
+        idx_uv.remove(axis)
+    else:  # 2D:
+        idx_uv.remove(2)
+    # Find default values:
+    if center is None:
+        center = tuple([d/2 for d in dim])
+    if radius is None:
+        radius = np.max((dim[idx_uv[0]], dim[idx_uv[1]])) / 4
+    if height is None and len(dim) == 3:  # only used for 3D!
+        height = dim[axis] / 2
+        assert height > 0, 'Height has to be a positive scalar value!'
+    assert len(dim) == len(center), 'Parameters dim and center must have the same dimensions!'
+    assert radius > 0, 'Radius has to be a positive scalar value!'
+    coords = np.indices(dim) + 0.5  # 0.5 to get to pixel/voxel center!
+    coords = coords - np.asarray(center)[(slice(None),) + (None,)*len(dim)]  # [:, None, None, None]
+    data = np.where(np.hypot(coords[idx_uv[0]], coords[idx_uv[1]]) <= radius, 1, 0)
+    if len(dim) == 3:  # 3D: Implement bounds above and below the disc:
+        height_shape = np.zeros_like(data)
+        bounds = [slice(None)] * 3  # i.e.: [:, :, :]
+        start = int(np.floor(center[axis] - height/2))
+        stop = int(np.ceil(center[axis] + height/2))
+        bounds[axis] = slice(start, stop)  # replace with actual bounds along disc symmetry axis
+        height_shape[tuple(bounds)] = 1
+        data *= height_shape
+    return Field(data=data, scale=scale, vector=False)
+
+
+def create_shape_ellipse(dim, center=None, width=None, height=None, axis=0, scale=1.0):
+    """Creates a Field object with the shape of an ellipse in 2D or 3D.
+
+    Attributes
+    ----------
+    dim : tuple
+        The dimensions of the grid.
+    center : tuple, optional
+        The center of the ellipse in pixel coordinates.
+    width : tuple, optional
+        The two half axes of the ellipse in pixel coordinates.
+    height : float, optional
+        The height of the ellipse in pixel coordinates. Unused if only 2D.
+    axis : int, optional
+        The orientation of the ellipses orthogonal axis. Only used in 3D case with z-axis as default.
+    scale: tuple of float
+        Scaling along the dimensions of the underlying data.
+
+    """
+    _log.debug('Calling __init__')
+    assert len(dim) in (2, 3), 'Ellipse can only be build in 2 or 3 dimensions!'
+    # Find indices of the disc plane axes:
+    idx_uv = [0, 1, 2]
+    if len(dim) == 3:  # 3D:
+        idx_uv.remove(axis)
+    else:  # 2D:
+        idx_uv.remove(2)
+    # Find default values:
+    if center is None:
+        center = tuple([d/2 for d in dim])
+    if width is None:
+        dim_uv = (dim[idx_uv[0]], dim[idx_uv[1]])
+        width = (np.max(dim_uv)*1/3, np.min(dim_uv)*2/3)
+    if height is None and len(dim) == 3:  # only used for 3D!
+        height = dim[axis] / 2
+    assert len(dim) == len(center), 'Parameters dim and center must have the same dimensions!'
+    assert len(width) == 2, 'Width has to contain the two half axes!'
+    coords = np.indices(dim) + 0.5  # 0.5 to get to pixel/voxel center!
+    coords = coords - np.asarray(center)[(slice(None),) + (None,)*len(dim)]  # [:, None, None, None]
+    distance = np.hypot(coords[idx_uv[1]] / (width[1]/2), coords[idx_uv[0]] / (width[0]/2))
+    data = np.where(distance <= 1, 1, 0)
+    if len(dim) == 3:  # 3D: Implement bounds above and below the disc:
+        height_shape = np.zeros_like(data)
+        bounds = [slice(None)] * 3  # i.e.: [:, :, :]
+        start = int(np.floor(center[axis] - height/2))
+        stop = int(np.ceil(center[axis] + height/2))
+        bounds[axis] = slice(start, stop)  # replace with actual bounds along disc symmetry axis
+        height_shape[tuple(bounds)] = 1
+        data *= height_shape
+    return Field(data=data, scale=scale, vector=False)
+
+
+def create_shape_ellipsoid(dim, center=None, width=None, scale=1.0):
+    """Creates a Field object with the shape of an ellipsoid in arbitrary dimensions.
+
+    Attributes
+    ----------
+    dim : tuple
+        The dimensions of the grid.
+    center : tuple, optional
+        The center of the ellipsoid in pixel coordinates.
+    width : tuple, optional
+        The half axes of the ellipsoid in pixel coordinates.
+    scale: tuple of float
+        Scaling along the dimensions of the underlying data.
+
+    """
+    _log.debug('Calling __init__')
+    if center is None:
+        center = tuple([d/2 for d in dim])
+    if width is None:
+        width = tuple([d/2 for d in dim])
+    assert len(dim) == len(center) == len(width), 'Parameters dim, center and width must have the same dimensions!'
+    coords = np.indices(dim) + 0.5  # 0.5 to get to pixel/voxel center!
+    coords = coords - np.asarray(center)[(slice(None),) + (None,)*len(dim)]  # [:, None, None, None]
+    distance = np.sqrt([(coords[i] / (width[i]/2))**2 for i in range(len(dim))])
+    data = np.where(distance <= 1, 1, 0)
+    return Field(data=data, scale=scale, vector=False)
+
+
+def create_shape_sphere(dim, center=None, radius=None, scale=1.0):
+    """Creates a Field object with the shape of a sphere in arbitrary dimensions.
+
+    Attributes
+    ----------
+    dim : tuple
+        The dimensions of the grid.
+    center : tuple, optional
+        The center of the sphere in pixel coordinates.
+    width : tuple, optional
+        The half axes of the sphere in pixel coordinates.
+    scale: tuple of float
+        Scaling along the dimensions of the underlying data.
+
+    """
+    _log.debug('Calling __init__')
+    if center is None:
+        center = tuple([d/2 for d in dim])
+    if radius is None:
+        radius = np.max(dim) / 4
+    assert len(dim) == len(center), 'Parameters dim and center must have the same dimensions!'
+    assert radius > 0, 'Radius has to be a positive scalar value!'
+    coords = np.indices(dim) + 0.5  # 0.5 to get to pixel/voxel center!
+    coords = coords - np.asarray(center)[(slice(None),) + (None,)*len(dim)]  # [:, None, None, None]
+    distance = np.sqrt([coords[i]**2 for i in range(len(dim))])
+    data = np.where(distance <= radius, 1, 0)
+    return Field(data=data, scale=scale, vector=False)
+
+
+def create_shape_filament(dim, pos=None, axis=0, scale=1.0):
+    """Creates a Field object with the shape of a filament in arbitrary dimension.
+
+    Parameters
+    ----------
+    dim : tuple
+        The dimensions of the grid.
+    pos : tuple, optional
+        The position of the filament in pixel coordinates. Has to be a tuple of len(dim) - 1 and denotes the
+        index positions along all axes aside from the specified `axis` along which the filament is placed.
+    axis :  int, optional
+        The orientation of the filament axis. Defaults to the first axis.
+    scale: tuple of float
+        Scaling along the dimensions of the underlying data.
+
+    """
+    _log.debug('Calling __init__')
+    assert len(dim) > 1, 'Only usable for multidimensional Fields (at least 2 dimensions)!'
+    if pos is None:
+        pos = tuple([d/2 for d in dim])
+    assert len(dim) == len(pos)-1, 'Position has to specify one less coordinates than dimensions!'
+    data = np.zeros(dim)
+    index = list(pos)
+    index.insert(axis, slice(None))  # [:] for the filament axis! pos everywhere else!
+    data[tuple(index)] = 1
+    return Field(data=data, scale=scale, vector=False)
+
+
+def create_shape_voxel(dim, pos=None, scale=1.0):
+    """Creates a Field object with the shape of a single voxel in arbitrary dimension.
+
+    Parameters
+    ----------
+    dim : tuple
+        The dimensions of the grid.
+    pos : tuple, optional
+        The position of the voxel.
+    scale: tuple of float
+        Scaling along the dimensions of the underlying data.
+
+    """
+    _log.debug('Calling __init__')
+    assert len(dim) > 1, 'Only usable for multidimensional Fields (at least 2 dimensions)!'
+    if pos is None:
+        pos = tuple([d/2 for d in dim])
+    assert len(dim) == len(pos), 'Parameters dim and pos must have the same dimensions!'
+    data = np.zeros(dim)
+    data[pos] = 1
+    return Field(data=data, scale=scale, vector=False)

src/empyre/fields/vectors.py

+# -*- coding: utf-8 -*-
+# Copyright 2016 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+
+
+import logging
+from numbers import Number
+
+import numpy as np
+
+from .field import Field
+
+
+__all__ = ['create_vector_homog', 'create_vector_vortex', 'create_vector_skyrmion', 'create_vector_singularity']
+
+_log = logging.getLogger(__name__)
+
+
+def create_vector_homog(dim, phi=0, theta=None, scale=1.0):
+    """Field subclass implementing a homogeneous vector field with 2 or 3 components in 2 or 3 dimensions.
+
+    Attributes
+    ----------
+    dim : tuple
+        The dimensions of the grid.
+    phi : float
+        The azimuthal angle. The default is 0, i.e., the vectors point in x direction.
+    theta : float, optional
+        The polar angle. If None (default), only two components will be created (corresponds to pi/2 in 3D, i.e., the
+        vectors are in the xy plane).
+    scale: tuple of float
+        Scaling along the dimensions of the underlying data.
+
+    """
+    _log.debug('Calling __init__')
+    assert len(dim) in (2, 3), 'Disc can only be build in 2 or 3 dimensions!'
+    assert isinstance(phi, Number), 'phi has to be an angle in radians!'
+    assert isinstance(theta, Number) or theta is None, 'theta has to be an angle in radians or None!'
+    if len(dim) == 2 and theta is None:  # 2D and 3rd component not explicitely wanted:
+        y_comp = np.ones(dim) * np.sin(phi)
+        x_comp = np.ones(dim) * np.cos(phi)
+        data = np.stack([x_comp, y_comp], axis=-1)
+    else:  # 3D, always have a third component:
+        if theta is None:
+            theta = np.pi/2  # xy-plane if not specified otherwise!
+        z_comp = np.ones(dim) * np.cos(theta)
+        y_comp = np.ones(dim) * np.sin(theta) * np.sin(phi)
+        x_comp = np.ones(dim) * np.sin(theta) * np.cos(phi)
+        data = np.stack([x_comp, y_comp, z_comp], axis=-1)
+    return Field(data=data, scale=scale, vector=True)
+
+
+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
+    ----------
+    dim : tuple
+        The dimensions of the grid.
+    center : tuple (N=2 or N=3), optional
+        The vortex center, given in 2D `(v, u)` or 3D `(z, y, x)`, where the perpendicular axis is discarded
+        (determined by the `axis` parameter). Is set to the center of the field of view if not specified.
+        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.
+    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 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. If this is `None`, only two components will be created for 2-dimensional vortices!
+    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 small.
+    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')
+    assert len(dim) in (2, 3), 'Vortex can only be build in 2 or 3 dimensions!'
+    # Find indices of the vortex plane axes:
+    idx_uv = [0, 1, 2]
+    if len(dim) == 3:  # 3D:
+        idx_uv.remove(axis)
+    else:  # 2D:
+        idx_uv.remove(2)
+    # 2D dimensions:
+    dim_uv = tuple([dim[i] for i in idx_uv])
+    # Find default values:
+    if center is None:
+        center = tuple([dim[i] / 2 for i in idx_uv])
+    elif len(center) == 3:  # if a 3D-center is given, just take the relevant coordinates:
+        center = list(center)
+        del center[axis]
+        center = tuple(center)
+    # 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]) - 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:
+        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)
+        u_comp = np.expand_dims(u_comp, axis=axis)
+        reps = [1, 1, 1]  # repetitions for tiling
+        reps[axis] = dim[axis]  # repeat along chosen axis
+        w_comp = np.tile(w_comp, reps)
+        v_comp = np.tile(v_comp, reps)
+        u_comp = np.tile(u_comp, reps)
+    if axis == 0:  # z-axis
+        z_comp = w_comp
+        y_comp = -v_comp
+        x_comp = -u_comp
+    elif axis == 1:  # y-axis
+        z_comp = v_comp
+        y_comp = w_comp
+        x_comp = u_comp
+    elif axis == 2:  # x-axis
+        z_comp = -v_comp
+        y_comp = -u_comp
+        x_comp = w_comp
+    else:
+        raise ValueError(f'{axis} is not a valid argument for axis (has to be 0, 1 or 2)')
+    data = np.stack([x_comp, y_comp, z_comp], axis=-1)
+    return Field(data=data, scale=scale, vector=True)
+
+
+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
+    ----------
+    dim : tuple
+        The dimensions of the grid.
+    center : tuple (N=2 or N=3), optional
+        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 center has to be between two pixels.
+    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 :  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
+    -------
+    field : `~.Field` object
+        The resulting vector field.
+
+    Notes
+    -----
+        To avoid singularities, the source center should lie between the pixel centers (which
+        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).
+
+        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 is described in a paper by Romming et al
+        (see DOI: 10.1103/PhysRevLett.114.177203s)
+
+    """
+
+    def _theta(r):
+        theta_1 = + np.arcsin(np.tanh((r + skyrm_d/2)/(wall_d/2)))
+        theta_2 = - np.arcsin(np.tanh((r - skyrm_d/2)/(wall_d/2)))
+        theta = theta_1 + theta_2
+        theta /= np.abs(theta).max() / np.pi
+        return theta
+
+    _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:
+        idx_uv.remove(axis)
+    else:  # 2D:
+        idx_uv.remove(2)
+    # 2D dimensions:
+    dim_uv = tuple([dim[i] for i in idx_uv])
+    # Find default values:
+    if skyrm_d is None:
+        skyrm_d = np.min(dim_uv) / 2
+    if wall_d is None:
+        wall_d = skyrm_d / 4
+    if center is None:
+        center = tuple([dim[i] / 2 for i in idx_uv])
+    elif len(center) == 3:  # if a 3D-center is given, just take the relevant coordinates:
+        center = list(center)
+        del center[axis]
+        center = tuple(center)
+    # Create skyrmion 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!
+    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 = 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:
+    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)
+        u_comp = np.expand_dims(u_comp, axis=axis)
+        reps = [1, 1, 1]  # repetitions for tiling
+        reps[axis] = dim[axis]  # repeat along chosen axis
+        w_comp = np.tile(w_comp, reps)
+        v_comp = np.tile(v_comp, reps)
+        u_comp = np.tile(u_comp, reps)
+    if axis == 0:  # z-axis
+        z_comp = w_comp
+        y_comp = -v_comp
+        x_comp = -u_comp
+    elif axis == 1:  # y-axis
+        z_comp = v_comp
+        y_comp = w_comp
+        x_comp = u_comp
+    elif axis == 2:  # x-axis
+        z_comp = -v_comp
+        y_comp = -u_comp
+        x_comp = w_comp
+    else:
+        raise ValueError(f'{axis} is not a valid argument for axis (has to be 0, 1 or 2)')
+    data = np.stack([x_comp, y_comp, z_comp], axis=-1)
+    return Field(data=data, scale=scale, vector=True)
+
+
+def create_vector_singularity(dim, center=None, scale=1.0):
+    """Create a 3-dimensional magnetic distribution of a homogeneous magnetized object.
+
+    Parameters
+    ----------
+    dim : tuple
+        The dimensions of the grid.
+    center : tuple (N=2 or N=3), optional
+        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.
+    scale: tuple of float
+        Scaling along the dimensions of the underlying data.
+
+    Returns
+    -------
+    field : `~.Field` object
+        The resulting vector field.
+
+    Notes
+    -----
+        To avoid singularities, the source center should lie between the pixel centers (which
+        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).
+
+        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:
+        center = tuple([d / 2 for d in dim])
+    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!
+    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))
+    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)

src/empyre/io/__init__.py

+# -*- coding: utf-8 -*-
+# Copyright 2020 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""Subpackage containing EMPyRe IO functionality for several EMPyRe classes."""
+
+
+from .io_field import *
+
+
+__all__ = []
+__all__.extend(io_field.__all__)
+
+
+del io_field

src/empyre/io/field_plugins/__init__.py

+# -*- coding: utf-8 -*-
+# Copyright 2020 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""Subpackage containing EMPyRe IO functionality for the Field class."""
+
+
+from . import llg, numpy, ovf, tec, text, vtk
+
+
+plugin_list = [llg, numpy, ovf, tec, text, vtk]
+__all__ = ['plugin_list']

src/empyre/io/field_plugins/llg.py

+# -*- coding: utf-8 -*-
+# Copyright 2020 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""IO plugin for LLG format."""
+
+
+import logging
+
+import numpy as np
+
+from ...fields.field import Field
+
+
+_log = logging.getLogger(__name__)
+
+file_extensions = ('.llg',)  # Recognised file extensions
+
+
+def reader(filename, scale=None, vector=None, **kwargs):
+    _log.debug('Call reader')
+    if vector is None:
+        vector = True
+    assert vector is True, 'Only vector fields can be read from the llg file format!'
+    SCALE = 1.0E-9 / 1.0E-2  # From cm to nm
+    data_columns = np.genfromtxt(filename, skip_header=2)
+    dim = tuple(np.genfromtxt(filename, dtype=int, skip_header=1, skip_footer=len(data_columns[:, 0])))
+    if scale is None:  # Otherwise overwrite!
+        stride_x = 1  # x varies fastest
+        stride_y = dim[2]  # one y step for dim[2] x steps
+        stride_z = dim[1] * dim[2]  # one z step for one x-y layer (dim[1]*dim[2])
+        scale_x = (data_columns[stride_x, 0] - data_columns[0, 0]) / SCALE  # first column varies in x
+        scale_y = (data_columns[stride_y, 1] - data_columns[0, 1]) / SCALE  # second column varies in y
+        scale_z = (data_columns[stride_z, 2] - data_columns[0, 2]) / SCALE  # third column varies in z
+        scale = (scale_z, scale_y, scale_x)
+    x_mag, y_mag, z_mag = data_columns[:, 3:6].T
+    data = np.stack((x_mag.reshape(dim), y_mag.reshape(dim), z_mag.reshape(dim)), axis=-1)
+    return Field(data, scale, vector=True)
+
+
+def writer(filename, field, **kwargs):
+    _log.debug('Call writer')
+    assert field.vector and len(field.dim) == 3, 'Only 3D vector fields can be saved to the llg file format!'
+    SCALE = 1.0E-9 / 1.0E-2  # from nm to cm
+    # Create 3D meshgrid and reshape it and the field into a list where x varies first:
+    zzz, yyy, xxx = (np.indices(field.dim) + 0.5) * np.reshape(field.scale, (3, 1, 1, 1)) * SCALE  # broadcast shape!
+    z_coord, y_coord, x_coord = np.ravel(zzz), np.ravel(yyy), np.ravel(xxx)  # Turn into vectors!
+    x_comp, y_comp, z_comp = field.comp  # Extract scalar field components!
+    x_vec, y_vec, z_vec = np.ravel(x_comp.data), np.ravel(y_comp.data), np.ravel(z_comp.data)  # Turn into vectors!
+    data = np.array([x_coord, y_coord, z_coord, x_vec, y_vec, z_vec]).T
+    # Save data to file:
+    with open(filename, 'w') as mag_file:
+        mag_file.write('LLGFileCreator: EMPyRe vector Field\n')
+        mag_file.write('    {:d}    {:d}    {:d}\n'.format(*field.dim))
+        mag_file.writelines('\n'.join('   '.join('{:7.6e}'.format(cell) for cell in row) for row in data))

src/empyre/io/field_plugins/numpy.py

+# -*- coding: utf-8 -*-
+# Copyright 2020 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""IO plugin for the numpy format."""
+
+
+import logging
+
+import numpy as np
+
+from ...fields.field import Field
+
+
+_log = logging.getLogger(__name__)
+
+file_extensions = ('.npy', '.npz')  # Recognised file extensions
+
+
+def reader(filename, scale=None, vector=None, **kwargs):
+    _log.debug('Call reader')
+    if vector is None:
+        vector = False
+    if scale is None:
+        scale = 1.0
+    return Field(np.load(filename, **kwargs), scale, vector)
+
+
+def writer(filename, field, **kwargs):
+    _log.debug('Call writer')
+    np.save(filename, field.data, **kwargs)

src/empyre/io/field_plugins/ovf.py

+# -*- coding: utf-8 -*-
+# Copyright 2020 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""IO plugin for simple text format."""
+
+
+import logging
+import os
+
+import numpy as np
+
+from ...fields.field import Field
+
+
+_log = logging.getLogger(__name__)
+
+file_extensions = ('.ovf', '.omf', '.ohf', 'obf')  # Recognised file extensions
+
+
+def reader(filename, scale=None, vector=None, segment=None, **kwargs):
+    """More info at:
+
+    http://math.nist.gov/oommf/doc/userguide11b2/userguide/vectorfieldformat.html
+    http://math.nist.gov/oommf/doc/userguide12a5/userguide/OVF_2.0_format.html
+
+    """
+    _log.debug('Calling reader')
+    if vector is None:
+        vector = True
+    assert vector is True, 'Only vector fields can be loaded from ovf-files!'
+    with open(filename, 'rb') as load_file:
+        line = load_file.readline()
+        assert line.startswith(b'# OOMMF')  # File has OVF format!
+        read_header, read_data = False, False
+        header = {'version': line.split()[-1].decode('utf-8')}
+        x_mag, y_mag, z_mag = [], [], []
+        data_mode = None
+        while True:
+            # --- READ START OF FILE OR IN BETWEEN SEGMENTS LINE BY LINE ---------------------------
+            if not read_header and not read_data:  # Start of file or between segments!
+                line = load_file.readline()
+                if line == b'':
+                    break  # End of file is reached!
+                if line.startswith(b'# Segment count'):
+                    seg_count = int(line.split()[-1])  # Total number of segments (often just 1)!
+                    seg_curr = 0  # Current segment (0: not in first segment, yet!)
+                    if seg_count > 1:  # If multiple segments, check if "segment" was set correctly:
+                        assert segment is not None, (f'Multiple ({seg_count}) segments were found! '
+                                                     'Chose one via the segment parameter!')
+                    elif segment is None:  # Only one segment AND parameter not set:
+                        segment = 1  # Default to the first/only segment!
+                    assert 0 < segment <= seg_count, (f'parameter segment={segment} out of bounds, '
+                                                      f'Use value between 1 and {seg_count}!')
+                    header['segment_count'] = seg_count
+                elif line.startswith(b'# Begin: Segment'):  # Segment start!
+                    seg_curr += 1
+                    if seg_curr > segment:
+                        break  # Stop reading the file!
+                elif line.startswith(b'# Begin: Header'):  # Header start!
+                    read_header = True
+                elif line.startswith(b'# Begin: Data'):  # Data start!
+                    read_data = True
+                    data_mode = ' '.join(line.decode('utf-8').split()[3:])
+                    assert data_mode in ['text', 'Text', 'Binary 4', 'Binary 8'], \
+                        f'Data mode {data_mode} is currently not supported by this reader!'
+                    assert header.get('meshtype') == 'rectangular', \
+                        'Only rectangular grids can be currently read!'
+            # --- READ HEADER LINE BY LINE ---------------------------------------------------------
+            elif read_header:
+                line = load_file.readline()
+                if line.startswith(b'# End: Header'):  # Header is done:
+                    read_header = False
+                    continue
+                line = line.decode('utf-8')  # Decode to use strings here!
+                line_list = line.split()
+                if '##' in line_list:  # Strip trailing comments:
+                    del line_list[line_list.index('##'):]
+                if len(line_list) <= 1:  # Just '#' or empty line:
+                    continue
+                key, value = line_list[1].strip(':'), ' '.join(line_list[2:])
+                if key not in header:  # Add new key, value pair if not existant:
+                    header[key] = value
+                elif key == 'Desc':  # Description can go over several lines:
+                    header['Desc'] = ' '.join([header['Desc'], value])
+            # --- READ DATA LINE BY LINE -----------------------------------------------------------
+            elif read_data:  # Currently in a data block:
+                if data_mode in ['text', 'Text']:  # Read data as text, line by line:
+                    line = load_file.readline()
+                    if line.startswith(b'# End: Data'):
+                        read_data = False  # Stop reading data and search for new segments!
+                        continue
+                    elif seg_curr < segment:  # Do nothing with the line if wrong segment!
+                        continue
+                    else:
+                        x, y, z = [float(i) for i in line.split()]
+                        x_mag.append(x)
+                        y_mag.append(y)
+                        z_mag.append(z)
+                elif 'Binary' in data_mode:  # Read data as binary, all bytes at the same time:
+                    # Currently every segment is read until the wanted one is processed. Only that one is returned!
+                    count = int(data_mode.split()[-1])  # Either 4 or 8!
+                    if header['version'] == '1.0':  # Big endian float:
+                        dtype = f'>f{count}'
+                    elif header['version'] == '2.0':  # Little endian float:
+                        dtype = f'<f{count}'
+                    test = np.fromfile(load_file, dtype=dtype, count=1)  # Read test byte!
+                    if count == 4:  # Binary 4:
+                        assert test == 1234567.0, 'Wrong test bytes!'
+                    elif count == 8:  # Binary 8:
+                        assert test == 123456789012345.0, 'Wrong test bytes!'
+                    dim = (int(header['znodes']), int(header['ynodes']), int(header['xnodes']))
+                    data_raw = np.fromfile(load_file, dtype=dtype, count=3*np.prod(dim))
+                    x_mag, y_mag, z_mag = data_raw[0::3], data_raw[1::3], data_raw[2::3]
+                    read_data = False  # Stop reading data and search for new segments (if any).
+        # --- READING DONE -------------------------------------------------------------------------
+        # Format after reading:
+        dim = (int(header['znodes']), int(header['ynodes']), int(header['xnodes']))
+        x_mag = np.asarray(x_mag).reshape(dim)
+        y_mag = np.asarray(y_mag).reshape(dim)
+        z_mag = np.asarray(z_mag).reshape(dim)
+        data = np.stack((x_mag, y_mag, z_mag), axis=-1) * float(header.get('valuemultiplier', 1))
+        if scale is None:
+            unit = header.get('meshunit', 'nm')
+            if unit == 'unspecified':
+                unit = 'nm'
+            _log.info(f'unit: {unit}')
+            unit_scale = {'m': 1e9, 'mm': 1e6, 'µm': 1e3, 'nm': 1}[unit]
+            xstep = float(header.get('xstepsize')) * unit_scale
+            ystep = float(header.get('ystepsize')) * unit_scale
+            zstep = float(header.get('zstepsize')) * unit_scale
+            scale = (zstep, ystep, xstep)
+    return Field(data, scale=scale, vector=True)
+
+
+def writer(filename, field, **kwargs):
+    _log.debug('Call writer')
+    assert field.vector and len(field.dim) == 3, 'Only 3D vector fields can be saved to ovf files!'
+    with open(filename, 'w') as save_file:
+        save_file.write('# OOMMF OVF 2.0\n')
+        save_file.write('# Segment count: 1\n')
+        save_file.write('# Begin: Segment\n')
+        # Write Header:
+        save_file.write('# Begin: Header\n')
+        name = os.path.split(filename)[1]
+        save_file.write(f'# Title: PYRAMID-VECTORDATA {name}\n')
+        save_file.write('# meshtype: rectangular\n')
+        save_file.write('# meshunit: nm\n')
+        save_file.write('# valueunit: A/m\n')
+        save_file.write('# valuemultiplier: 1.\n')
+        save_file.write('# xmin: 0.\n')
+        save_file.write('# ymin: 0.\n')
+        save_file.write('# zmin: 0.\n')
+        save_file.write(f'# xmax: {field.scale[2] * field.dim[2]}\n')
+        save_file.write(f'# ymax: {field.scale[1] * field.dim[1]}\n')
+        save_file.write(f'# zmax: {field.scale[0] * field.dim[0]}\n')
+        save_file.write('# xbase: 0.\n')
+        save_file.write('# ybase: 0.\n')
+        save_file.write('# zbase: 0.\n')
+        save_file.write(f'# xstepsize: {field.scale[2]}\n')
+        save_file.write(f'# ystepsize: {field.scale[1]}\n')
+        save_file.write(f'# zstepsize: {field.scale[0]}\n')
+        save_file.write(f'# xnodes: {field.dim[2]}\n')
+        save_file.write(f'# ynodes: {field.dim[1]}\n')
+        save_file.write(f'# znodes: {field.dim[0]}\n')
+        save_file.write('# End: Header\n')
+        # Write data:
+        save_file.write('# Begin: Data Text\n')
+        x_mag, y_mag, z_mag = field.comp
+        x_mag = x_mag.data.ravel()
+        y_mag = y_mag.data.ravel()
+        z_mag = z_mag.data.ravel()
+        for i in range(np.prod(field.dim)):
+            save_file.write(f'{x_mag[i]:g} {y_mag[i]:g} {z_mag[i]:g}\n')
+        save_file.write('# End: Data Text\n')
+        save_file.write('# End: Segment\n')

src/empyre/io/field_plugins/tec.py

+# -*- coding: utf-8 -*-
+# Copyright 2020 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""IO plugin for simple text format."""
+
+
+import logging
+import re
+
+import numpy as np
+
+from ...fields.field import Field
+from ...utils.misc import interp_to_regular_grid
+
+
+_log = logging.getLogger(__name__)
+
+file_extensions = ('.tec',)  # Recognised file extensions
+
+
+def reader(filename, scale=None, vector=None, **kwargs):
+    assert isinstance(scale, tuple), 'The scale must be a tuple, each entry corresponding to one grid dimensions!'
+    _log.debug('Call reader')
+    if vector is None:
+        vector = True
+    assert vector is True, 'Only vector fields can be loaded from tec-files!'
+    with open(filename, 'r') as mag_file:
+        lines = mag_file.readlines()  # Read in lines!
+        match = re.search(R'N=(\d+)', lines[2])  # Extract number of points from third line!
+        if match:
+            n_points = int(match.group(1))
+        else:
+            raise IOError('File does not seem to match .tec format!')
+    n_head, n_foot = 3, len(lines) - (3 + n_points)
+    # Read in data:
+    data_raw = np.genfromtxt(filename, skip_header=n_head, skip_footer=n_foot)
+    if scale is None:
+        raise ValueError('For the interpolation of unstructured grids, the `scale` parameter is required!')
+    data = interp_to_regular_grid(data_raw[:, :3], data_raw[:, 3:], scale, **kwargs)
+    return Field(data, scale, vector=vector)
+
+
+def writer(filename, field, **kwargs):
+    _log.debug('Call writer')
+    raise NotImplementedError('A writer for this extension is not yet implemented!')

src/empyre/io/field_plugins/text.py

+# -*- coding: utf-8 -*-
+# Copyright 2020 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""IO plugin for simple text format."""
+
+
+import logging
+import ast
+
+import numpy as np
+
+from ...fields.field import Field
+
+
+_log = logging.getLogger(__name__)
+
+file_extensions = ('.txt',)  # Recognised file extensions
+
+
+def reader(filename, scale=None, vector=None, **kwargs):
+    _log.debug('Call reader')
+    if vector is None:
+        vector = False
+    assert vector is False, 'Only scalar 2D fields can currently be read with this file reader!'
+    with open(filename, 'r') as load_file:  # Read data:
+        empyre_format = load_file.readline().startswith('EMPYRE-FORMAT')
+        if empyre_format:  # File has EMPyRe structure:
+            scale = ast.literal_eval(load_file.readline()[8:-4])  # [8:-4] takes just the scale string!
+            data = np.loadtxt(filename, delimiter='\t', skiprows=2)  # skips header!
+        else:  # Try default with provided kwargs:
+            scale = 1.0 if scale is None else scale  # Set default if not provided!
+            data = np.loadtxt(filename, **kwargs)
+    return Field(data, scale=scale, vector=False)
+
+
+def writer(filename, field, with_header=True, **kwargs):
+    _log.debug('Call writer')
+    assert not field.vector, 'Vector fields can currently not be saved to text!'
+    assert len(field.dim) == 2, 'Only 2D fields can currenty be saved to text!'
+    if with_header:  # write header:
+        with open(filename, 'w') as save_file:
+            save_file.write('EMPYRE-FORMAT\n')
+            save_file.write(f'scale = {field.scale} nm\n')
+        save_kwargs = {'fmt': '%7.6e', 'delimiter': '\t'}
+    else:
+        save_kwargs = kwargs
+    with open(filename, 'ba') as save_file:  # the 'a' is for append!
+        np.savetxt(save_file, field.data, **save_kwargs)

src/empyre/io/field_plugins/vtk.py

+# -*- coding: utf-8 -*-
+# Copyright 2020 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""IO plugin for simple text format."""
+
+
+import logging
+from numbers import Number
+
+import numpy as np
+
+from ...fields.field import Field
+from ...utils.misc import interp_to_regular_grid, restrict_points
+from ...vis import colors
+
+
+_log = logging.getLogger(__name__)
+
+file_extensions = ('.vtk',)  # Recognised file extensions
+
+
+def reader(filename, scale=None, vector=None, bounds=None, **kwargs):
+    """More infos at:
+
+    overview: https://docs.enthought.com/mayavi/mayavi/data.html
+    writing: https://vtk.org/Wiki/VTK/Writing_VTK_files_using_python
+    format: https://www.vtk.org/wp-content/uploads/2015/04/file-formats.pdf
+
+    """
+    _log.debug('Calling reader')
+    try:
+        from tvtk.api import tvtk
+    except ImportError:
+        _log.error('This extension recquires the tvtk package!')
+        return
+    if vector is None:
+        vector = True
+    # Setting up reader:
+    reader = tvtk.DataSetReader(file_name=filename, read_all_scalars=True, read_all_vectors=vector)
+    reader.update()
+    # Getting output:
+    output = reader.output
+    assert output is not None, 'File reader could not find data or file "{}"!'.format(filename)
+    # Reading points and vectors:
+    if isinstance(output, tvtk.ImageData):  # tvtk.StructuredPoints is a subclass of tvtk.ImageData!
+        # Connectivity: implicit; described by: 3D data array and spacing along each axis!
+        _log.info('geometry: ImageData')
+        # Load relevant information from output (reverse to get typical Python order z,y,x):
+        dim = output.dimensions[::-1]
+        origin = output.origin[::-1]
+        spacing = output.spacing[::-1]
+        _log.info(f'dim: {dim}, origin: {origin}, spacing: {spacing}')
+        assert len(dim) == 3, 'Data has to be three-dimensional!'
+        if scale is None:
+            scale = tuple(spacing)
+        if vector:  # Extract vector compontents and create magnitude array:
+            vector_array = np.asarray(output.point_data.vectors, dtype=np.float)
+            x_mag, y_mag, z_mag = vector_array.T
+            data = np.stack((x_mag.reshape(dim), y_mag.reshape(dim), z_mag.reshape(dim)), axis=-1)
+        else:  # Extract scalar data and create magnitude array:
+            scalar_array = np.asarray(output.point_data.scalars, dtype=np.float)
+            data = scalar_array.reshape(dim)
+    elif isinstance(output, tvtk.RectilinearGrid):
+        # Connectivity: implicit; described by: 3D data array and 1D array of spacing for each axis!
+        _log.info('geometry: RectilinearGrid')
+        raise NotImplementedError('RectilinearGrid is currently not supported!')
+    elif isinstance(output, tvtk.StructuredGrid):
+        # Connectivity: implicit; described by: 3D data array and 3D position arrays for each axis!
+        _log.info('geometry: StructuredGrid')
+        raise NotImplementedError('StructuredGrid is currently not supported!')
+    elif isinstance(output, tvtk.PolyData):
+        # Connectivity: explicit; described by: x, y, z positions of vertices and arrays of surface cells!
+        _log.info('geometry: PolyData')
+        raise NotImplementedError('PolyData is currently not supported!')
+    elif isinstance(output, tvtk.UnstructuredGrid):
+        # Connectivity: explicit; described by: x, y, z positions of vertices and arrays of volume cells!
+        _log.info('geometry: UnstructuredGrid')
+        # Load relevant information from output:
+        point_array = np.asarray(output.points, dtype=np.float)
+        if vector:
+            data_array = np.asarray(output.point_data.vectors, dtype=np.float)
+        else:
+            data_array = np.asarray(output.point_data.scalars, dtype=np.float)
+        if scale is None:
+            raise ValueError('For the interpolation of unstructured grids, the `scale` parameter is required!')
+        elif isinstance(scale, Number):  # Scale is the same for each dimension x, y, z!
+            scale = (scale,) * 3
+        elif isinstance(scale, tuple):
+            assert len(scale) == 3, f'Each dimension (z, y, x) needs a scale, but {scale} was given!'
+        # Crop data to required range, if necessary
+        if bounds is not None:
+            _log.info('Restrict data')
+            point_array, data_array = restrict_points(point_array, data_array, bounds)
+        data = interp_to_regular_grid(point_array, data_array, scale, **kwargs)
+    else:
+        raise TypeError('Data type of {} not understood!'.format(output))
+    return Field(data, scale, vector=True)
+
+
+def writer(filename, field, **kwargs):
+    _log.debug('Call writer')
+    assert len(field.dim) == 3, 'Currently only 3D fields can be saved to vtk!'
+    try:
+        from tvtk.api import tvtk, write_data
+    except ImportError:
+        _log.error('This extension recquires the tvtk package!')
+        return
+    # Create dataset:
+    origin = (0, 0, 0)
+    spacing = (field.scale[2], field.scale[1], field.scale[0])
+    dimensions = (field.dim[2], field.dim[1], field.dim[0])
+    sp = tvtk.StructuredPoints(origin=origin, spacing=spacing, dimensions=dimensions)
+    # Fill with data from field:
+    if field.vector:  # Handle vector fields:
+        # Put vector components in corresponding array:
+        vectors = field.data.reshape(-1, 3)
+        sp.point_data.vectors = vectors
+        sp.point_data.vectors.name = 'vectors'
+        # 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()))
+        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'
+        point_colors.from_array(rgb)
+        sp.point_data.scalars = point_colors
+        sp.point_data.scalars.name = 'colors'
+    else:  # Handle scalar fields:
+        scalars = field.data.ravel()
+        sp.point_data.scalars = scalars
+        sp.point_data.scalars.name = 'scalars'
+    # Write the data to file:
+    write_data(sp, filename)

src/empyre/io/io_field.py

+# -*- coding: utf-8 -*-
+# Copyright 2020 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""IO functionality for Field objects."""
+
+
+import logging
+import os
+
+from ..fields.field import Field
+from .field_plugins import plugin_list
+
+
+__all__ = ['load_field', 'save_field']
+_log = logging.getLogger(__name__)
+
+
+def load_field(filename, scale=None, vector=None, **kwargs):
+    """Load supported file into a :class:`~..fields.Field` instance.
+
+    The function loads the file according to the extension:
+        SCALAR???
+        - hdf5 for HDF5.  # TODO: You can use comp_pos here!!!
+        - EMD Electron Microscopy Dataset format (also HDF5).
+        - npy or npz for numpy formats.
+
+        PHASEMAP
+        - hdf5 for HDF5.
+        - rpl for Ripple (useful to export to Digital Micrograph).
+        - dm3 and dm4 for Digital Micrograph files.
+        - unf for SEMPER unf binary format.
+        - txt format.
+        - npy or npz for numpy formats.
+        - Many image formats such as png, tiff, jpeg...
+
+        VECTOR
+        - hdf5 for HDF5.
+        - EMD Electron Microscopy Dataset format (also HDF5).
+        - llg format.
+        - ovf format.
+        - npy or npz for numpy formats.
+
+    Any extra keyword is passed to the corresponsing reader. For available options see their individual documentation.
+
+    Parameters
+    ----------
+    filename:  str
+        The filename to be loaded.
+    scale: tuple of float, optional
+        Scaling along the dimensions of the underlying data. Default is 1.
+    vector: bool, optional
+        True if the field should be a vector field, False if it should be interpreted as a scalar field (default).
+
+    Returns
+    -------
+    field: `Field`
+        A `Field` object containing the loaded data.
+
+    Notes
+    -----
+    Falls back to HyperSpy routines for loading data, make sure it is installed if you need the full capabilities.
+
+    """
+    _log.debug('Calling load_field')
+    extension = os.path.splitext(filename)[1]
+    for plugin in plugin_list:  # Iterate over all plugins:
+        if extension in plugin.file_extensions:  # Check if extension is recognised:
+            return plugin.reader(filename, scale=scale, vector=vector, **kwargs)
+    # If nothing was found, try HyperSpy
+    _log.debug('Using HyperSpy')
+    try:
+        import hyperspy.api as hs
+    except ImportError:
+        _log.error('This extension recquires the hyperspy package!')
+        return
+    comp_pos = kwargs.pop('comp_pos', -1)
+    return Field.from_signal(hs.load(filename, **kwargs), scale=scale, vector=vector, comp_pos=comp_pos)
+
+
+def save_field(filename, field, **kwargs):
+    """Saves the Field  in the specified format.
+
+    The function gets the format from the extension:
+        - hdf5 for HDF5.
+        - EMD Electron Microscopy Dataset format (also HDF5).
+        - npy or npz for numpy formats.
+
+    If no extension is provided, 'hdf5' is used. Most formats are saved with the HyperSpy package (internally the field
+    is first converted to a HyperSpy Signal.
+    Each format accepts a different set of parameters. For details see the specific format documentation.
+
+    Parameters
+    ----------
+    filename : str, optional
+        Name of the file which the Field is saved into. The extension determines the saving procedure.
+
+    """
+    """Saves the phasemap in the specified format.
+
+    The function gets the format from the extension:
+        - hdf5 for HDF5.
+        - rpl for Ripple (useful to export to Digital Micrograph).
+        - unf for SEMPER unf binary format.
+        - txt format.
+        - Many image formats such as png, tiff, jpeg...
+
+    If no extension is provided, 'hdf5' is used. Most formats are
+    saved with the HyperSpy package (internally the phasemap is first
+    converted to a HyperSpy Signal.
+
+    Each format accepts a different set of parameters. For details
+    see the specific format documentation.
+
+    Parameters
+    ----------
+    filename: str, optional
+        Name of the file which the phasemap is saved into. The extension
+        determines the saving procedure.
+    save_mask: boolean, optional
+        If True, the `mask` is saved, too. For all formats, except HDF5, a separate file will
+        be created. HDF5 always saves the `mask` in the metadata, independent of this flag. The
+        default is False.
+    save_conf: boolean, optional
+        If True, the `confidence` is saved, too. For all formats, except HDF5, a separate file
+        will be created. HDF5 always saves the `confidence` in the metadata, independent of
+        this flag. The default is False
+    pyramid_format: boolean, optional
+        Only used for saving to '.txt' files. If this is True, the grid spacing is saved
+        in an appropriate header. Otherwise just the phase is written with the
+        corresponding `kwargs`.
+
+    """
+    _log.debug('Calling save_field')
+    extension = os.path.splitext(filename)[1]
+    for plugin in plugin_list:  # Iterate over all plugins:
+        if extension in plugin.file_extensions:  # Check if extension is recognised:
+            plugin.writer(filename, field, **kwargs)
+            return
+    # If nothing was found, try HyperSpy:
+    _log.debug('Using HyperSpy')
+    field.to_signal().save(filename, **kwargs)

src/empyre/utils/__init__.py

+# -*- coding: utf-8 -*-
+# Copyright 2020 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""Subpackage containing utility functionality."""
+
+
+from .quaternion import *
+
+
+__all__ = []
+__all__.extend(quaternion.__all__)
+
+
+del quaternion

src/empyre/utils/misc.py

+# -*- coding: utf-8 -*-
+# Copyright 2020 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""This module provides the miscellaneous helper functions."""
+
+
+import logging
+import itertools
+from time import time
+
+import numpy as np
+from tqdm import tqdm
+from scipy.spatial import cKDTree, qhull
+from scipy.interpolate import LinearNDInterpolator
+
+
+__all__ = ['levi_civita', 'interp_to_regular_grid', 'restrict_points']
+_log = logging.getLogger(__name__)
+
+
+def levi_civita(i, j, k):
+    _log.debug('Calling levi_civita')
+    return (j-i) * (k-i) * (k-j) / (np.abs(j-i) * np.abs(k-i) * np.abs(k-j))
+
+
+def interp_to_regular_grid(points, values, scale, scale_factor=1, step=1, convex=True, distance_upper_bound=None):
+    """Interpolate values on points to regular grid.
+
+    Parameters
+    ----------
+    points : np.ndarray, (N, 3)
+        Array of points, describing the location of the values that should be interpolated. Three columns x, y, z!
+    values : np.ndarray, (N, c)
+        Array of values that should be interpolated to the new grid. `c` is the number of components (`1` for scalar
+        fields, `3` for normal 3D vector fields).
+    scale : tuple of 3 ints
+        Scale along each of the 3 spatial dimensions. Usually given in nm.
+    scale_factor : float, optional
+        Additional scaling factor that should be used if the original points are not described on a nm-scale. Use this
+        to convert from nm of `scale` to the unit of the points. By default 1.
+    step : int, optional
+        If this is bigger than 1 (the default), only every `step` point is taken into account. Can speed up calculation,
+        but you'll lose accuracy of your interpolation.
+    convex : bool, optional
+        By default True. If this is set to False, additional measures are taken to find holes in the point cloud.
+        WARNING: this is an experimental feature that should be used with caution!
+    distance_upper_bound: float, optional
+        Only used if `convex=False`. Set the upper bound, determining if a point of the new (interpolated) grid is too
+        far away from any original point. They are assumed to be in a "hole" and their values are set to zero. Set this
+        value in nm, it will be converted to the local unit of the original points internally. If not set and
+        `convex=True`, double of the the mean of `scale` is calculated and used (can be troublesome if the scales vary
+        drastically).
+
+    Returns
+    -------
+    interpolation: np.ndarray
+        Interpolated grid with shape `(zdim, ydim, xdim)` for scalar and `(zdim, ydim, xdim, ncomp)` for vector field.
+
+    """
+    _log.debug('Calling interpolate_to_regular_grid')
+    z_uniq = np.unique(points[:, 2])
+    _log.info(f'unique positions along z: {len(z_uniq)}')
+    #  Translate scale in local units (not necessarily nm), taken care of with `scale_factor`:
+    scale = tuple([s * scale_factor for s in scale])
+    # Determine the size of the point cloud of irregular coordinates:
+    x_min, x_max = points[:, 0].min(), points[:, 0].max()
+    y_min, y_max = points[:, 1].min(), points[:, 1].max()
+    z_min, z_max = points[:, 2].min(), points[:, 2].max()
+    x_diff, y_diff, z_diff = np.ptp(points[:, 0]), np.ptp(points[:, 1]), np.ptp(points[:, 2])
+    _log.info(f'x-range: {x_min:.2g} <-> {x_max:.2g} ({x_diff:.2g})')
+    _log.info(f'y-range: {y_min:.2g} <-> {y_max:.2g} ({y_diff:.2g})')
+    _log.info(f'z-range: {z_min:.2g} <-> {z_max:.2g} ({z_diff:.2g})')
+    # Determine dimensions from given grid spacing a:
+    dim = tuple(np.round(np.asarray((z_diff/scale[0], y_diff/scale[1], x_diff/scale[2]), dtype=int)))
+    assert all(dim) > 0, f'All dimensions of dim={dim} need to be > 0, please adjust the scale accordingly!'
+    z = z_min + scale[0] * (np.arange(dim[0]) + 0.5)  # +0.5: shift to pixel center!
+    y = y_min + scale[1] * (np.arange(dim[1]) + 0.5)  # +0.5: shift to pixel center!
+    x = x_min + scale[2] * (np.arange(dim[2]) + 0.5)  # +0.5: shift to pixel center!
+    # Create points for new Euclidian grid; fliplr for (x, y, z) order:
+    points_euc = np.fliplr(np.asarray(list(itertools.product(z, y, x))))
+    # Make values 2D (if not already); double .T so that a new axis is added at the END (n, 1):
+    values = np.atleast_2d(values.T).T
+    # Prepare interpolated grid:
+    interpolation = np.empty(dim+(values.shape[-1],), dtype=np.float)
+    _log.info(f'Dimensions of new grid: {interpolation.shape}')
+    # Calculate the Delaunay triangulation (same for every component of multidim./vector fields):
+    _log.info('Start Delaunay triangulation...')
+    tick = time()
+    triangulation = qhull.Delaunay(points[::step])
+    tock = time()
+    _log.info(f'Delaunay triangulation complete (took {tock-tick:.2f} s)!')
+    # Perform the interpolation for each column of `values`:
+    for i in tqdm(range(values.shape[-1])):
+        # Create interpolator for the given triangulation and the values of the current column:
+        interpolator = LinearNDInterpolator(triangulation, values[::step, i], fill_value=0)
+        # Interpolate:
+        interpolation[..., i] = interpolator(points_euc).reshape(dim)
+    # If NOT convex, we have to check for additional holes in the structure (EXPERIMENTAL):
+    if not convex:  # Only necessary if the user expects holes in the (-> nonconvex) distribution:
+        # Create k-dimensional tree for queries:
+        _log.info('Create cKDTree...')
+        tick = time()
+        tree = cKDTree(points)
+        tock = time()
+        _log.info(f'cKDTree creation complete (took {tock-tick:.2f} s)!')
+        # Query the tree for nearest neighbors, x: points to query, k: number of neighbors, p: norm
+        # to use (here: 2 - Euclidean), distance_upper_bound: maximum distance that is searched!
+        if distance_upper_bound is None:  # Take the mean of the scale for the upper bound:
+            distance_upper_bound = 2 * np.mean(scale)  # NOTE: could be problematic for wildly varying scale numbers.
+        else:  # Convert to local scale:
+            distance_upper_bound *= scale_factor
+        _log.info('Start cKDTree neighbour query...')
+        tick = time()
+        data, leafsize = tree.query(x=points_euc, k=1, p=2, distance_upper_bound=distance_upper_bound)
+        tock = time()
+        _log.info(f'cKDTree neighbour query complete (took {tock-tick:.2f} s)!')
+        # Create boolean mask that determines which interpolation points have no neighbor near enough:
+        mask = np.isinf(data).reshape(dim)  # Points further away than upper bound were marked 'inf'!
+        _log.info(f'{np.sum(mask)} of {points_euc.shape[0]} points were assumed to be in holes of the point cloud!')
+        # Set these points to zero (NOTE: This can take a looooong time):
+        interpolation[mask, :] = 0
+    return np.squeeze(interpolation)
+
+
+def restrict_points(point_array, data_array, bounds):
+    """Restrict range of point_array and data_array
+
+    Parameters
+    ----------
+    points_array : np.ndarray, (N, 3)
+        Array of points, describing the location of the values that should be interpolated. Three columns x, y, z!
+    data_array : np.ndarray, (N, c)
+        Array of values that should be interpolated to the new grid. `c` is the number of components (`1` for scalar
+        fields, `3` for normal 3D vector fields).
+    bounds : tuple of 6 values
+        Restrict data range to given bounds, x0, x1, y0, y1, z0, z1.
+
+    Returns
+    -------
+    point_restricted: np.ndarray
+        Cut out of the array of points inside the bounds, describing the location of the values that should be
+        interpolated. Three columns x, y, z!
+    value_restricted: np.ndarray
+        Cut out of the array of values inside the bounds, describing the location of the values that should be
+        interpolated. Three columns x, y, z!
+
+    """
+    point_restricted = []
+    data_restricted = []
+    for i, pos in enumerate(point_array):
+        if bounds[0] <= pos[0] <= bounds[1]:
+            if bounds[2] <= pos[1] <= bounds[3]:
+                if bounds[4] <= pos[2] <= bounds[5]:
+                    point_restricted.append(pos)
+                    data_restricted.append(data_array[i])
+    point_restricted = np.array(point_restricted)
+    data_restricted = np.array(data_restricted)
+    return point_restricted, data_restricted

pyramid/quaternion.py → src/empyre/utils/quaternion.py

-# -*- coding: utf-8 -*-
-# Copyright 2014 by Forschungszentrum Juelich GmbH
-# Author: J. Caron
-#
-"""This module provides the :class:`~.Quaternion` class which can be used for rotations."""
-
-import logging
-
-import numpy as np
-
-__all__ = ['Quaternion']
-
-
-class Quaternion(object):
-    """Class representing a rotation expressed by a quaternion.
-
-    A quaternion is a four-dimensional description of a rotation which can also be described by
-    a rotation vector (`v1`, `v2`, `v3`) and a rotation angle :math:`\theta`. The four components
-    are calculated to:
-
-    .. math::
-
-        w = \cos(\theta/2)
-        x = v_1 \cdot \sin(\theta/2)
-        y = v_2 \cdot \sin(\theta/2)
-        z = v_3 \cdot \sin(\theta/2)
-
-    Use the :func:`~.from_axisangle` and :func:`~.to_axisangle` to convert to axis-angle
-    representation and vice versa. Quaternions can be multiplied by other quaternions, which
-    results in a new rotation or with a vector, which results in a rotated vector.
-
-    Attributes
-    ----------
-    values : float
-        The four quaternion values `w`, `x`, `y`, `z`.
-
-    """
-
-    NORM_TOLERANCE = 1E-6
-
-    _log = logging.getLogger(__name__ + '.Quaternion')
-
-    @property
-    def conj(self):
-        """The conjugate of the quaternion, representing a tilt in opposite direction."""
-        w, x, y, z = self.values
-        return Quaternion((w, -x, -y, -z))
-
-    @property
-    def matrix(self):
-        """The rotation matrix representation of the quaternion."""
-        w, x, y, z = self.values
-        return np.array([[1 - 2 * (y ** 2 + z ** 2), 2 * (x * y - w * z), 2 * (x * z + w * y)],
-                         [2 * (x * y + w * z), 1 - 2 * (x ** 2 + z ** 2), 2 * (y * z - w * x)],
-                         [2 * (x * z - w * y), 2 * (y * z + w * x), 1 - 2 * (x ** 2 + y ** 2)]])
-
-    def __init__(self, values):
-        self._log.debug('Calling __init__')
-        self.values = values
-        self._normalize()
-        self._log.debug('Created ' + str(self))
-
-    def __mul__(self, other):  # self * other
-        self._log.debug('Calling __mul__')
-        if isinstance(other, Quaternion):  # Quaternion multiplication
-            return self.dot_quat(self, other)
-        elif len(other) == 3:  # vector multiplication
-            q_vec = Quaternion((0,) + tuple(other))
-            q = self.dot_quat(self.dot_quat(self, q_vec), self.conj)
-            return q.values[1:]
-
-    def dot_quat(self, q1, q2):
-        """Multiply two :class:`~.Quaternion` objects to create a new one (always normalized).
-
-        Parameters
-        ----------
-        q1, q2 : :class:`~.Quaternion`
-            The quaternion which should be multiplied.
-
-        Returns
-        -------
-        quaternion : :class:`~.Quaternion`
-            The resulting quaternion.
-
-        """
-        self._log.debug('Calling dot_quat')
-        w1, x1, y1, z1 = q1.values
-        w2, x2, y2, z2 = q2.values
-        w = w1 * w2 - x1 * x2 - y1 * y2 - z1 * z2
-        x = w1 * x2 + x1 * w2 + y1 * z2 - z1 * y2
-        y = w1 * y2 + y1 * w2 + z1 * x2 - x1 * z2
-        z = w1 * z2 + z1 * w2 + x1 * y2 - y1 * x2
-        return Quaternion((w, x, y, z))
-
-    def _normalize(self):
-        self._log.debug('Calling _normalize')
-        mag2 = np.sum(n ** 2 for n in self.values)
-        if abs(mag2 - 1.0) > self.NORM_TOLERANCE:
-            mag = np.sqrt(mag2)
-            self.values = tuple(n / mag for n in self.values)
-
-    @classmethod
-    def from_axisangle(cls, vector, theta):
-        """Create a quaternion from an axis-angle representation
-
-        Parameters
-        ----------
-        vector : :class:`~numpy.ndarray` (N=3)
-            Vector around which the rotation is executed.
-        theta : float
-            Rotation angle.
-
-        Returns
-        -------
-        quaternion : :class:`~.Quaternion`
-            The resulting quaternion.
-
-        """
-        cls._log.debug('Calling from_axisangle')
-        x, y, z = vector
-        theta /= 2.
-        w = np.cos(theta)
-        x *= np.sin(theta)
-        y *= np.sin(theta)
-        z *= np.sin(theta)
-        return cls((w, x, y, z))
-
-    def to_axisangle(self):
-        """Convert the quaternion to axis-angle-representation.
-
-        Returns
-        -------
-        vector, theta : :class:`~numpy.ndarray` (N=3), float
-            Vector around which the rotation is executed and rotation angle.
-
-        """
-        self._log.debug('Calling to_axisangle')
-        w, x, y, z = self.values
-        theta = 2.0 * np.arccos(w)
-        return np.array((x, y, z)), theta
+# -*- coding: utf-8 -*-
+# Copyright 2020 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""This module provides the :class:`~.Quaternion` class which can be used for rotations."""
+
+
+import logging
+
+import numpy as np
+
+
+__all__ = ['Quaternion']
+
+
+class Quaternion(object):
+    R"""Class representing a rotation expressed by a quaternion.
+
+    A quaternion is a four-dimensional description of a rotation which can also be described by
+    a rotation vector (`v1`, `v2`, `v3`) and a rotation angle :math:`\theta`. The four components
+    are calculated to:
+
+    .. math::
+
+        w = \cos(\theta/2)
+        x = v_1 \cdot \sin(\theta/2)
+        y = v_2 \cdot \sin(\theta/2)
+        z = v_3 \cdot \sin(\theta/2)
+
+    Use the :func:`~.from_axisangle` and :func:`~.to_axisangle` to convert to axis-angle
+    representation and vice versa. Quaternions can be multiplied by other quaternions, which
+    results in a new rotation or with a vector, which results in a rotated vector.
+
+    Attributes
+    ----------
+    values : float
+        The four quaternion values `w`, `x`, `y`, `z`.
+
+    """
+
+    NORM_TOLERANCE = 1E-6
+
+    _log = logging.getLogger(__name__ + '.Quaternion')
+
+    @property
+    def conj(self):
+        """The conjugate of the quaternion, representing a tilt in opposite direction."""
+        w, x, y, z = self.values
+        return Quaternion((w, -x, -y, -z))
+
+    @property
+    def matrix(self):
+        """The rotation matrix representation of the quaternion."""
+        w, x, y, z = self.values
+        return np.array([[1 - 2 * (y ** 2 + z ** 2), 2 * (x * y - w * z), 2 * (x * z + w * y)],
+                         [2 * (x * y + w * z), 1 - 2 * (x ** 2 + z ** 2), 2 * (y * z - w * x)],
+                         [2 * (x * z - w * y), 2 * (y * z + w * x), 1 - 2 * (x ** 2 + y ** 2)]])
+
+    def __init__(self, values):
+        self._log.debug('Calling __init__')
+        self.values = values
+        self._normalize()
+        self._log.debug('Created ' + str(self))
+
+    def __mul__(self, other):  # self * other
+        self._log.debug('Calling __mul__')
+        if isinstance(other, Quaternion):  # Quaternion multiplication
+            return self.dot_quat(self, other)
+        elif len(other) == 3:  # vector multiplication (Caution: normalises!)
+            q_vec = Quaternion((0,) + tuple(other))
+            q = self.dot_quat(self.dot_quat(self, q_vec), self.conj)
+            return q.values[1:]
+
+    def _normalize(self):
+        self._log.debug('Calling _normalize')
+        mag2 = np.sum([n ** 2 for n in self.values])
+        if abs(mag2 - 1.0) > self.NORM_TOLERANCE:
+            mag = np.sqrt(mag2)
+            self.values = tuple(n / mag for n in self.values)
+
+    def dot_quat(self, q1, q2):
+        """Multiply two :class:`~.Quaternion` objects to create a new one (always normalized).
+
+        Parameters
+        ----------
+        q1, q2 : :class:`~.Quaternion`
+            The quaternion which should be multiplied.
+
+        Returns
+        -------
+        quaternion : :class:`~.Quaternion`
+            The resulting quaternion.
+
+        """
+        self._log.debug('Calling dot_quat')
+        w1, x1, y1, z1 = q1.values
+        w2, x2, y2, z2 = q2.values
+        w = w1 * w2 - x1 * x2 - y1 * y2 - z1 * z2
+        x = w1 * x2 + x1 * w2 + y1 * z2 - z1 * y2
+        y = w1 * y2 + y1 * w2 + z1 * x2 - x1 * z2
+        z = w1 * z2 + z1 * w2 + x1 * y2 - y1 * x2
+        return Quaternion((w, x, y, z))
+
+    @classmethod
+    def from_axisangle(cls, vector, theta):
+        """Create a quaternion from an axis-angle representation
+
+        Parameters
+        ----------
+        vector : :class:`~numpy.ndarray` (N=3)
+            Vector around which the rotation is executed.
+        theta : float
+            Rotation angle.
+
+        Returns
+        -------
+        quaternion : :class:`~.Quaternion`
+            The resulting quaternion.
+
+        """
+        cls._log.debug('Calling from_axisangle')
+        x, y, z = vector
+        theta /= 2.
+        w = np.cos(theta)
+        x *= np.sin(theta)
+        y *= np.sin(theta)
+        z *= np.sin(theta)
+        return cls((w, x, y, z))
+
+    def to_axisangle(self):
+        """Convert the quaternion to axis-angle-representation.
+
+        Returns
+        -------
+        vector, theta : :class:`~numpy.ndarray` (N=3), float
+            Vector around which the rotation is executed and rotation angle.
+
+        """
+        self._log.debug('Calling to_axisangle')
+        w, x, y, z = self.values
+        theta = 2.0 * np.arccos(w)
+        return np.array((x, y, z)), theta