environment.yml

-# Create environment with "conda env create", update with "conda env update"!
-# This file contains the top level dependencies of the empyre project!
-# requirements.txt also contains sub-dependencies, generated by "pip freeze > requirements.txt"!
-# Add new dependencies here, then "conda env update", then "pip freeze > requirements.txt"!
-# To see if compatible upgrades are available for the current packages, use "conda upgrade --all"!
-# When packages are deprecated/deleted, it may be best to recreate the environment from scratch!
-
-name: empyre
-
-channels:
-  - defaults  # Default conda channels, on top to keep highest priority!
-  - conda-forge  # Used for hyperspy, pyFFTW!
-
-dependencies:
-  # Basic:
-  - python=3.7
-  - setuptools
-  - numpy=1.17
-  - scipy=1.3
-  - tqdm=4.36
-  - scikit-image=0.15
-  # File IO:
-  - hyperspy=1.5
-  - hyperspy-gui-ipywidgets=1.2
-  #- h5py=2.9  # TODO: dependency of hyperspy? Not needed here?
-  # Plotting and colors:
-  - matplotlib=3.1
-  - Pillow=6.1
-  - cmocean=2.0
-  # 3D plotting:
-  - mayavi=4.7
-  - ipyevents=0.7
-  - ipywidgets=7.5
-  # Testing:
-  - pytest=5.0
-  - pytest-cov=2.7
-  - pytest-flake8=1.0
-  - pytest-runner=5.1
-  - pyroma  # Checks setup.py for cheese!
-  #- pytest-mpl=0.10  # Needed for testing hyperspy! # TODO: Use for pyramid/plotting library, too!
-  - coverage=4.5
-  # Documentation:
-  - sphinx=2.4
-  - numpydoc=0.9
-  - sphinx_rtd_theme=0.4
-  # IPython and notebooks:
-  - ipython=7.7
-  - jupyter=1.0
-  - nb_conda=2.2
-  #- ptvsd=4.3  # Cell debugging in VS Code
-  - rope=0.16  # for refactoring in VS Code
-  # TODO: - ipywidgets
-  # TODO: Add back GUI dependencies!
-  # TODO: Get Jutil from gitlab (currently doesn't work, git and cygwin don't play nice,...
-  # TODO: ...because one is Unix, ond is Windows).
-  # Fast computation:
-  - pyFFTW=0.11
-  # TODO: ? - pathos  # pathos.multiprocessing uses dill instead of pickle
-  # PIP installations:
-  #- pip=19.0
-#  - pip:
-#    # ALSO NEEDS JUTIL!
-#    - "git+ssh://gitlab@iffgit.fz-juelich.de/unger/jutil.git"
-  # Misc.:
-  - nodejs=12.8  # TODO: Needs to be fixed to prevent errors! Check again if needed in future!

pyproject.toml

+[build-system]
+requires = ["hatchling"]
+build-backend = "hatchling.build"
+
+[project]
+name = "empyre"
+dynamic = ["version"]
+description = "Electron Microscopy Python Reconstruction"
+readme = "README.rst"
+license = "GPL-3.0-or-later"
+requires-python = ">=3.6"
+authors = [
+    { name = "Jan Caron", email = "j.caron@fz-juelich.de" },
+]
+keywords = [
+    "Electron Microscopy",
+    "Inverse Problem Solving",
+    "Model-based Reconstrution",
+]
+classifiers = [
+    "Development Status :: 3 - Alpha",
+    "Intended Audience :: Developers",
+    "Intended Audience :: Science/Research",
+    "License :: OSI Approved :: GNU General Public License v3 (GPLv3)",
+    "Operating System :: OS Independent",
+    "Programming Language :: Python :: 3.6",
+    "Programming Language :: Python :: 3.7",
+    "Programming Language :: Python :: 3.8",
+    "Programming Language :: Python :: 3.9",
+    "Programming Language :: Python :: 3.10",
+    "Topic :: Scientific/Engineering",
+]
+dependencies = [
+    "matplotlib >= 3.2",
+    "numpy >= 1.17",
+    "Pillow",
+    "scikit-image",
+    "scipy",
+    "tqdm",
+]
+
+[project.optional-dependencies]
+3d = [
+    "mayavi >= 4.7",
+]
+all = [
+    "cmocean",
+    "pyfftw",
+]
+colors = [
+    "cmocean",
+]
+docs = [
+    "nbsphinx",
+    "nbsphinx_link",
+    "numpydoc",
+    "pandoc",
+    "sphinx",
+    "sphinx_rtd_theme",
+]
+fftw = [
+    "pyfftw",
+]
+io = [
+    "hyperspy",
+    "tvtk",
+]
+tests = [
+    "coverage",
+    "pyroma",
+    "pytest",
+    "pytest-cov",
+    "pytest-flake8",
+    "pytest-runner",
+]
+
+[project.urls]
+Homepage = "https://iffgit.fz-juelich.de/empyre/empyre"
+
+[tool.hatch.version]
+path = "src/empyre/__init__.py"
+
+[tool.hatch.build.targets.sdist]
+include = [
+    "/src",
+    "/tests",
+]

setup.cfg

@@ -2,75 +2,10 @@
 # setup.cfg
 # :copyright: Copyright 2020 Jan Caron
 
-
-# CONFIGURATION FOR SETUP.PY:
-
-[metadata]
-name = empyre
-version = 0.2.1
-author = Jan Caron
-author-email = j.caron@fz-juelich.de
-description = Electron Microscopy Python Reconstruction
-long-description = file: README.rst
-url = https://iffgit.fz-juelich.de/empyre/empyre
-license = GPLv3
-classifiers =
-    Development Status :: 3 - Alpha
-    Intended Audience :: Developers
-    Intended Audience :: Science/Research
-    License :: OSI Approved :: GNU General Public License v3 (GPLv3)
-    Programming Language :: Python :: 3.7
-    Topic :: Scientific/Engineering
-    Operating System :: OS Independent
-keywords =
-    Electron Microscopy
-    Inverse Problem Solving
-    Model-based Reconstrution
-
-[options]
-zip_safe = False
-include_package_data = True
-package_dir =
-    =src
-packages = find:
-python_requires = >=3.7
-setup_requires =
-    setuptools
-tests_require =
-    coverage
-    pytest
-    pytest-cov
-    pytest-flake8
-    pytest-runner
-install_requires =
-    numpy >= 1.17
-    matplotlib >= 3
-    scikit-image
-    tqdm
-    scipy
-    Pillow
-
-[options.packages.find]
-where = src
-
-[options.extras_require]
-io =
-    hyperspy
-    tvtk
-fftw =
-    pyfftw
-colors =
-    cmocean
-3d =
-    mayavi >= 4.7
-all =
-    tvtk
-    pyfftw
-    cmocean
-    mayavi >= 4.7
 # TODO: pip install .[io] is not working because hyperspy depends on trait which has no wheels on PyPI at the moment...
 # TODO: See https://github.com/hyperspy/hyperspy/issues/2315 and https://github.com/enthought/traits/issues/357
 # TODO: Add hyperspy back as soon (if?) this is resolved... Until then: install hyperspy with conda!
+# TODO: mayavi has the same exact problem and tvtk is apparently also a problem...
 
 # CONFIGURATION FOR TESTING:
 
@@ -86,13 +21,20 @@ omit =
 [flake8]
 max-line-length = 120
 ignore =
-    E402  # module import not at top of file
-    E124  # closing bracket does not match visual indentation
-    E125  # continuation line with same indent as next logical line
-    E226  # missing whitespace around arithmetic operator
-    W503  # line break before binary operator
-    W504  # line break after binary operator
-    E741  # do not use variables named ‘l’, ‘O’, or ‘I’
+    # module import not at top of file
+    E402
+    # closing bracket does not match visual indentation
+    E124
+    # continuation line with same indent as next logical line
+    E125
+    # missing whitespace around arithmetic operator
+    E226
+    # line break before binary operator
+    W503
+    # line break after binary operator
+    W504
+    # do not use variables named ‘l’, ‘O’, or ‘I’
+    E741
 per-file-ignores =
     */__init__.py: F401, F403, F405, F821
     # F401: module imported but unused

setup.py

-#!python
-# coding=utf-8
-"""Setup for testing, building, distributing and installing the 'EMPyRe'-package"""
-
-
-import os
-import subprocess
-from setuptools import setup, config
-
-
-# Read version from setup.py:
-version = config.read_configuration('setup.cfg')['metadata']['version']
-# Get current git revision:
-try:
-    git_rev = subprocess.check_output(['git', 'rev-parse', 'HEAD']).strip().decode()
-except Exception:
-    git_rev = "???"
-# Write both to version.py:
-print(R'writing src\empyre\version.py')
-with open(os.path.join(os.path.dirname(__file__), 'src', 'empyre', 'version.py'), 'w') as vfile:
-    vfile.write('# -*- coding: utf-8 -*-\n' +
-                '""""This file was automatically generated by `setup.py`"""\n' +
-                f'version = "{version}"\n' +
-                f'git_revision = "{git_rev}"\n')
-
-# Run setup (reads metadata & options from setup.cfg):
-print(R'running setup.py')
-setup()

src/empyre/__init__.py

@@ -11,16 +11,13 @@ from . import models
 from . import reconstruct
 from . import vis
 from . import utils
-from .version import version as __version__
-from .version import git_revision as __git_revision__
+
+__version__ = "0.3.4"
+__git_revision__ = "undefined"
 
 import logging
 _log = logging.getLogger(__name__)
-_log.info(f'Imported EMPyRe V-{__version__} GIT-{__git_revision__}')
+_log.info(f'Imported EMPyRe V-{__version__}')
 del logging
 
-
 __all__ = ['fields', 'io', 'models', 'reconstruct', 'vis', 'utils']
-
-
-del version

src/empyre/fields/shapes.py

@@ -20,7 +20,7 @@ _log = logging.getLogger(__name__)
 # TODO: TEST!!!
 
 
-def create_shape_slab(dim, center=None, width=None, scale=1):
+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
@@ -51,7 +51,7 @@ def create_shape_slab(dim, center=None, width=None, scale=1):
     return Field(data=data, scale=scale, vector=False)
 
 
-def create_shape_disc(dim, center=None, radius=None, height=None, axis=0, scale=1):
+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
@@ -102,7 +102,7 @@ def create_shape_disc(dim, center=None, radius=None, height=None, axis=0, scale=
     return Field(data=data, scale=scale, vector=False)
 
 
-def create_shape_ellipse(dim, center=None, width=None, height=None, axis=0, scale=1):
+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
@@ -154,7 +154,7 @@ def create_shape_ellipse(dim, center=None, width=None, height=None, axis=0, scal
     return Field(data=data, scale=scale, vector=False)
 
 
-def create_shape_ellipsoid(dim, center=None, width=None, scale=1):
+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
@@ -182,7 +182,7 @@ def create_shape_ellipsoid(dim, center=None, width=None, scale=1):
     return Field(data=data, scale=scale, vector=False)
 
 
-def create_shape_sphere(dim, center=None, radius=None, scale=1):
+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
@@ -211,7 +211,7 @@ def create_shape_sphere(dim, center=None, radius=None, scale=1):
     return Field(data=data, scale=scale, vector=False)
 
 
-def create_shape_filament(dim, pos=None, axis=0, scale=1):
+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
@@ -239,7 +239,7 @@ def create_shape_filament(dim, pos=None, axis=0, scale=1):
     return Field(data=data, scale=scale, vector=False)
 
 
-def create_shape_voxel(dim, pos=None, scale=1):
+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

src/empyre/fields/vectors.py

@@ -17,8 +17,8 @@ __all__ = ['create_vector_homog', 'create_vector_vortex', 'create_vector_skyrmio
 _log = logging.getLogger(__name__)
 
 
-def create_vector_homog(dim, phi=0, theta=None, scale=1):
-    """Field subclass implementing a homogeneous vector field with 3 components in 2 or 3 dimensions.
+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
     ----------
@@ -37,11 +37,13 @@ def create_vector_homog(dim, phi=0, theta=None, scale=1):
     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 theta is 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:
+    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)
@@ -49,9 +51,7 @@ def create_vector_homog(dim, phi=0, theta=None, scale=1):
     return Field(data=data, scale=scale, vector=True)
 
 
-def create_vector_vortex(dim, center=None, core_r=0, oop_r=None, axis=0, scale=1):
-    # TODO: oop and core_r documentation!!! General description!
-    # TODO: CHIRALITY
+def create_vector_vortex(dim, center=None, phi_0=np.pi/2, oop_r=None, oop_sign=1, core_r=0, axis=0, scale=1.0):
     """Field subclass implementing a vortex vector field with 3 components in 2 or 3 dimensions.
 
     Attributes
@@ -64,7 +64,31 @@ def create_vector_vortex(dim, center=None, core_r=0, oop_r=None, axis=0, scale=1
         The vortex center should be between two pixels to avoid singularities.
     axis :  int, optional
         The orientation of the vortex axis. The default is 0 and corresponds to the z-axis. Is ignored if dim is
-        only 2D. To invert chirality, multiply the resulting Field object by -1.  # TODO: WRONG?
+        only 2D.
+    phi_0 : float, optional
+        Angular offset that allows all arrows to be rotated simultaneously. The default value is `pi/2` and corresponds
+        to an anti-clockwise rotation, while a value of `-pi/2` corresponds to a clockwise rotation. `0` and `pi` lead
+        to arrows that all point in/outward.
+    oop_r : float, optional
+        Radius of a potential out-of-plane ("oop") component in the vortex center. If `None`, the vortex is purely
+        in-plane, which is the default. Set this to a number (in pixels) that determines the 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')
@@ -87,15 +111,19 @@ def create_vector_vortex(dim, center=None, core_r=0, oop_r=None, axis=0, scale=1
     # Create vortex plane (2D):
     coords_uv = np.indices(dim_uv) + 0.5  # 0.5 to get to pixel/voxel center!
     coords_uv = coords_uv - np.asarray(center, dtype=float)[:, None, None]  # Shift by center!
-    phi = np.arctan2(coords_uv[0], coords_uv[1]) - np.pi / 2
+    phi = np.arctan2(coords_uv[0], coords_uv[1]) - phi_0
     rr = np.hypot(coords_uv[0], coords_uv[1])
     rr_clip = np.clip(rr - core_r, a_min=0, a_max=None)  # negative inside core_r (clipped to 0), positive outside!
+    # Check if a core should be constructed (oop_r != 0):
     if oop_r is None:
         w_comp = np.zeros(dim_uv)
     else:
-        w_comp = 1 - 2/np.pi * np.arcsin(np.tanh(np.pi*rr_clip/oop_r))  # orthogonal: 1 inside, towards 0 outside!
-    v_comp = np.ones(dim_uv) * np.sin(phi) * np.sqrt(1 - w_comp)
-    u_comp = np.ones(dim_uv) * np.cos(phi) * np.sqrt(1 - w_comp)
+        assert oop_r >= 0, 'oop_r has to be a positive number!'
+        assert oop_sign in (1, -1), 'Sign of the out-of-plane component has to be either +1 or -1!'
+        w_comp = 1 - 2/np.pi * np.arcsin(np.tanh(np.pi*rr_clip/(oop_r+1E-30)))  # orthogonal: 1 inside, to 0 outside!
+        w_comp *= oop_sign * w_comp
+    v_comp = np.ones(dim_uv) * np.sin(phi) * np.sqrt(1 - np.abs(w_comp))
+    u_comp = np.ones(dim_uv) * np.cos(phi) * np.sqrt(1 - np.abs(w_comp))
     if len(dim) == 3:  # Expand to 3D:
         w_comp = np.expand_dims(w_comp, axis=axis)
         v_comp = np.expand_dims(v_comp, axis=axis)
@@ -123,7 +151,7 @@ def create_vector_vortex(dim, center=None, core_r=0, oop_r=None, axis=0, scale=1
     return Field(data=data, scale=scale, vector=True)
 
 
-def create_vector_skyrmion(dim, center=None, phi_0=0, skyrm_d=None, wall_d=None, axis=0, scale=1):
+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
@@ -137,20 +165,24 @@ def create_vector_skyrmion(dim, center=None, phi_0=0, skyrm_d=None, wall_d=None,
     phi_0 : float, optional
         Angular offset switching between Neel type (0 [default] or pi) or Bloch type (+/- pi/2)
         skyrmions.
+    core_sign : {1, -1}, optional
+        Can be `1` (default) or `-1` and determines if the skyrmion core is aligned parallel (`+1`) or antiparallel
+        (`-1`) to the chosen symmetry axis.
     skyrm_d : float, optional
         Diameter of the skyrmion. Defaults to half of the smaller dimension perpendicular to the
         skyrmion axis if not specified.
     wall_d : float, optional
         Diameter of the domain wall of the skyrmion. Defaults to `skyrm_d / 4` if not specified.
-    axis :  {'z', '-z', 'y', '-y', 'x', '-x'}, optional # TODO: NUMBERS!! See vortex
-        The orientation of the skyrmion axis. The default is 'z'. Negative values invert skyrmion
-        core direction.
+    axis :  int, optional
+        The orientation of the vortex axis. The default is 0 and corresponds to the z-axis. Is ignored if dim is
+        only 2D.
+    scale: tuple of float
+        Scaling along the dimensions of the underlying data.
 
     Returns
     -------
-    amplitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3)
-        The magnetic distribution as a tuple of the 3 components in
-        `x`-, `y`- and `z`-direction on the 3-dimensional grid.
+    field : `~.Field` object
+        The resulting vector field.
 
     Notes
     -----
@@ -159,9 +191,8 @@ def create_vector_skyrmion(dim, center=None, phi_0=0, skyrm_d=None, wall_d=None,
         coordinates (e.g. coordinate 1 lies between the first and the second pixel).
 
         Skyrmion wall width is dependant on exchange stiffness  A [J/m] and anisotropy K [J/m³]
-        The out-of-plane magnetization at the domain wall can be described as:
-        Mz = -Ms * tanh(x/w)  # TODO: Instead ROMER paper
-        w = sqrt(A/K)
+        The out-of-plane magnetization at the domain wall is described in a paper by Romming et al
+        (see DOI: 10.1103/PhysRevLett.114.177203s)
 
     """
 
@@ -174,6 +205,7 @@ def create_vector_skyrmion(dim, center=None, phi_0=0, skyrm_d=None, wall_d=None,
 
     _log.debug('Calling create_vector_skyrmion')
     assert len(dim) in (2, 3), 'Skyrmion can only be build in 2 or 3 dimensions!'
+    assert core_sign in (1, -1), 'Sign of the out-of-plane component has to be either +1 or -1!'
     # Find indices of the skyrmion plane axes:
     idx_uv = [0, 1, 2]
     if len(dim) == 3:  # 3D:
@@ -199,7 +231,7 @@ def create_vector_skyrmion(dim, center=None, phi_0=0, skyrm_d=None, wall_d=None,
     rr = np.hypot(coords_uv[0], coords_uv[1])
     phi = np.arctan2(coords_uv[0], coords_uv[1]) - phi_0
     theta = _theta(rr)
-    w_comp = np.cos(theta)
+    w_comp = core_sign * np.cos(theta)
     v_comp = np.sin(theta) * np.sin(phi)
     u_comp = np.sin(theta) * np.cos(phi)
     # Expansion to 3D if necessary and component shuffling:
@@ -230,7 +262,7 @@ def create_vector_skyrmion(dim, center=None, phi_0=0, skyrm_d=None, wall_d=None,
     return Field(data=data, scale=scale, vector=True)
 
 
-def create_vector_singularity(dim, center=None, scale=1):
+def create_vector_singularity(dim, center=None, scale=1.0):
     """Create a 3-dimensional magnetic distribution of a homogeneous magnetized object.
 
     Parameters
@@ -241,16 +273,13 @@ def create_vector_singularity(dim, center=None, scale=1):
         The source center, given in 2D `(v, u)` or 3D `(z, y, x)`, where the perpendicular axis
         is discarded. Is set to the center of the field of view if not specified.
         The source center has to be between two pixels.
-    axis :  {'z', '-z', 'y', '-y', 'x', '-x'}, optional  # TODO: NUMBERS!
-        The orientation of the source axis. The default is 'z'. Negative values invert the source
-        to a sink.  # TODO: wat...?
-    # TODO: scale!
+    scale: tuple of float
+        Scaling along the dimensions of the underlying data.
 
     Returns
     -------
-    amplitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3)
-        The magnetic distribution as a tuple of the 3 components in
-        `x`-, `y`- and `z`-direction on the 3-dimensional grid.
+    field : `~.Field` object
+        The resulting vector field.
 
     Notes
     -----
@@ -258,7 +287,10 @@ def create_vector_singularity(dim, center=None, scale=1):
         reside at coordinates with _.5 at the end), i.e. integer values should be used as center
         coordinates (e.g. coordinate 1 lies between the first and the second pixel).
 
-    """  # TODO: What does negating do here? Senke / Quelle (YES IT DOES!)? swell and sink? ISSUE! INVITE PEOPLE!!!
+        Per default, all arrows point outwards, negate the resulting `Field` object with a minus after creation to let
+        them point inwards.
+
+    """
     _log.debug('Calling create_vector_singularity')
     # Find default values:
     if center is None:
@@ -266,10 +298,12 @@ def create_vector_singularity(dim, center=None, scale=1):
     assert len(dim) == len(center), f"Length of dim ({len(dim)}) and center ({len(center)}) don't match!"
     # Setup coordinates, shape is (c, z, y, x), if 3D, or (c, y, x), if 2D (c: components):
     coords = np.indices(dim) + 0.5  # 0.5 to get to pixel/voxel center!
-    center = np.asarray(center, dtype=float)
-    bc_shape = (len(dim,),) + (1,)*len(dim)  # Shape for broadcasting, (1,1,1,3) for 3D, (1,1,2) for 2D!
-    coords = coords - center.reshape(bc_shape)  # Shift by center (append 1s for broadcasting)!
+    bc_shape = (len(dim,),) + (1,)*len(dim)  # Shape for broadcasting, (3,1,1,1) for 3D, (2,1,1) for 2D!
+    coords = coords - np.reshape(center, bc_shape)  # Shift by center (append 1s for broadcasting)!
     rr = np.sqrt(np.sum([coords[i]**2 for i in range(len(dim))], axis=0))
-    data = coords / (rr + 1E-30)  # Normalise amplitude (keep direction), rr (z,y,x) is broadcasted to data (c,z,y,x)!
-    data = data.T  # (c,z,y,x) -> (x,y,z,c)
+    coords /= rr + 1E-30  # Normalise amplitude (keep direction), rr (z,y,x) is broadcasted to data (c,z,y,x)!
+    # coords has components listed in wrong order (z, y, x) in the first dimension, we need (x, y, z) in the last:
+    coords = np.flip(coords, axis=0)
+    # Finally, the data has the coordinate axis at the end, not at the front:
+    data = np.moveaxis(coords, 0, -1)  # (c,z,y,x) -> (z,y,x,c)
     return Field(data=data, scale=scale, vector=True)

src/empyre/io/field_plugins/numpy.py

@@ -22,7 +22,7 @@ def reader(filename, scale=None, vector=None, **kwargs):
     if vector is None:
         vector = False
     if scale is None:
-        scale = 1
+        scale = 1.0
     return Field(np.load(filename, **kwargs), scale, vector)
 
 

src/empyre/io/field_plugins/ovf.py

@@ -18,7 +18,7 @@ _log = logging.getLogger(__name__)
 file_extensions = ('.ovf', '.omf', '.ohf', 'obf')  # Recognised file extensions
 
 
-def reader(filename, scale=None, vector=True, segment=None, **kwargs):
+def reader(filename, scale=None, vector=None, segment=None, **kwargs):
     """More info at:
 
     http://math.nist.gov/oommf/doc/userguide11b2/userguide/vectorfieldformat.html
@@ -26,6 +26,8 @@ def reader(filename, scale=None, vector=True, segment=None, **kwargs):
 
     """
     _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()
@@ -140,7 +142,7 @@ def writer(filename, field, **kwargs):
         save_file.write('# Begin: Segment\n')
         # Write Header:
         save_file.write('# Begin: Header\n')
-        name = os.path.split(os.path.split(filename)[1])
+        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')

src/empyre/io/field_plugins/tec.py

@@ -7,7 +7,6 @@
 
 import logging
 import re
-from numbers import Number
 
 import numpy as np
 
@@ -20,10 +19,11 @@ _log = logging.getLogger(__name__)
 file_extensions = ('.tec',)  # Recognised file extensions
 
 
-def reader(filename, dim, scale=None, vector=None, **kwargs):
+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 = False
+        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!
@@ -37,12 +37,8 @@ def reader(filename, dim, scale=None, vector=None, **kwargs):
     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!')
-    elif isinstance(scale, Number):  # Scale is the same for each dimension!
-        scale = (scale,) * len(dim)
-    elif isinstance(scale, tuple):
-        assert len(scale) == len(dim), f'Each of the {len(dim)} dimensions needs a scale, but {scale} was given!'
     data = interp_to_regular_grid(data_raw[:, :3], data_raw[:, 3:], scale, **kwargs)
-    return Field(data, scale, vector=False)
+    return Field(data, scale, vector=vector)
 
 
 def writer(filename, field, **kwargs):

src/empyre/io/field_plugins/text.py

@@ -29,7 +29,7 @@ def reader(filename, scale=None, vector=None, **kwargs):
             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 if scale is None else scale  # Set default if not provided!
+            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)
 

src/empyre/io/field_plugins/vtk.py

@@ -11,7 +11,7 @@ from numbers import Number
 import numpy as np
 
 from ...fields.field import Field
-from ...utils.misc import interp_to_regular_grid
+from ...utils.misc import interp_to_regular_grid, restrict_points
 from ...vis import colors
 
 
@@ -20,7 +20,7 @@ _log = logging.getLogger(__name__)
 file_extensions = ('.vtk',)  # Recognised file extensions
 
 
-def reader(filename, scale=None, vector=None, **kwargs):
+def reader(filename, scale=None, vector=None, bounds=None, **kwargs):
     """More infos at:
 
     overview: https://docs.enthought.com/mayavi/mayavi/data.html
@@ -43,8 +43,6 @@ def reader(filename, scale=None, vector=None, **kwargs):
     output = reader.output
     assert output is not None, 'File reader could not find data or file "{}"!'.format(filename)
     # Reading points and vectors:
-    tvtk.ImageData
-
     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')
@@ -86,10 +84,14 @@ def reader(filename, scale=None, vector=None, **kwargs):
             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!
-            scale = (scale,) * len(dim)
+        elif isinstance(scale, Number):  # Scale is the same for each dimension x, y, z!
+            scale = (scale,) * 3
         elif isinstance(scale, tuple):
-            assert len(scale) == len(dim), f'Each of the {len(dim)} dimensions needs a scale, but {scale} was given!'
+            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))
@@ -106,7 +108,7 @@ def writer(filename, field, **kwargs):
         return
     # Create dataset:
     origin = (0, 0, 0)
-    spacing = (field.scale[2], field.dim[1], field.scale[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:
@@ -118,7 +120,10 @@ def writer(filename, field, **kwargs):
         # Calculate colors:
         x_mag, y_mag, z_mag = field.comp
         magvec = np.asarray((x_mag.data.ravel(), y_mag.data.ravel(), z_mag.data.ravel()))
-        rgb = colors.CMAP_CIRCULAR_DEFAULT.rgb_from_vector(magvec)
+        cmap = kwargs.pop('cmap', None)
+        if cmap is None:
+            cmap = colors.cmaps.cyclic_cubehelix
+        rgb = cmap.rgb_from_vector(magvec)
         point_colors = tvtk.UnsignedIntArray()
         point_colors.number_of_components = 3
         point_colors.name = 'colors'

src/empyre/io/io_field.py

@@ -21,7 +21,7 @@ def load_field(filename, scale=None, vector=None, **kwargs):
 
     The function loads the file according to the extension:
         SCALAR???
-        - hdf5 for HDF5.
+        - hdf5 for HDF5.  # TODO: You can use comp_pos here!!!
         - EMD Electron Microscopy Dataset format (also HDF5).
         - npy or npz for numpy formats.
 
@@ -64,17 +64,18 @@ def load_field(filename, scale=None, vector=None, **kwargs):
     """
     _log.debug('Calling load_field')
     extension = os.path.splitext(filename)[1]
-    for i, plugin in enumerate(plugin_list):  # Iterate over all plugins:
+    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, vector, **kwargs)
-    # If nothing was found, try HyperSpy:
+            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
-    return Field.from_signal(hs.load(filename, **kwargs), scale, vector)
+    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):
@@ -132,7 +133,7 @@ def save_field(filename, field, **kwargs):
     """
     _log.debug('Calling save_field')
     extension = os.path.splitext(filename)[1]
-    for i, plugin in enumerate(plugin_list):  # Iterate over all plugins:
+    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

src/empyre/utils/misc.py

@@ -15,7 +15,7 @@ from scipy.spatial import cKDTree, qhull
 from scipy.interpolate import LinearNDInterpolator
 
 
-__all__ = ['levi_civita', 'interp_to_regular_grid']
+__all__ = ['levi_civita', 'interp_to_regular_grid', 'restrict_points']
 _log = logging.getLogger(__name__)
 
 
@@ -24,7 +24,7 @@ def levi_civita(i, j, k):
     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):
+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
@@ -45,6 +45,12 @@ def interp_to_regular_grid(points, values, scale, scale_factor=1, step=1, convex
     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
     -------
@@ -66,10 +72,11 @@ def interp_to_regular_grid(points, values, scale, scale_factor=1, step=1, convex
     _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, y_diff/scale, x_diff/scale), dtype=int)))
-    x = x_min + scale * (np.arange(dim[2]) + 0.5)  # +0.5: shift to pixel center!
-    y = y_min + scale * (np.arange(dim[1]) + 0.5)  # +0.5: shift to pixel center!
-    z = z_min + scale * (np.arange(dim[0]) + 0.5)  # +0.5: shift to pixel center!
+    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):
@@ -92,12 +99,61 @@ def interp_to_regular_grid(points, values, scale, scale_factor=1, step=1, convex
     # 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!
-        data, leafsize = tree.query(x=points_euc, k=1, p=2, distance_upper_bound=2*np.mean(scale))
+        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'!
-        for i in tqdm(range(values.shape[-1])):  # Set these points to zero (NOTE: This can take a looooong time):
-            interpolation[..., i].ravel()[mask.ravel()] = 0
+        _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

src/empyre/utils/quaternion.py

@@ -73,7 +73,7 @@ class Quaternion(object):
 
     def _normalize(self):
         self._log.debug('Calling _normalize')
-        mag2 = np.sum(n ** 2 for n in self.values)
+        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)

src/empyre/vis/colors.py

@@ -5,8 +5,7 @@
 
 """This module provides a number of custom colormaps, which also have capabilities for 3D plotting.
 If this is the case, the :class:`~.Colormap3D` colormap class is a parent class. In `cmaps`, a
-number of specialised colormaps is available for convenience. If the default for circular colormaps
-(used for 3D plotting) should be changed, set it via `CMAP_CIRCULAR_DEFAULT`.
+number of specialised colormaps is available for convenience.
 For general questions about colors see:
 http://www.poynton.com/PDFs/GammaFAQ.pdf
 http://www.poynton.com/PDFs/ColorFAQ.pdf
@@ -28,7 +27,7 @@ from .tools import use_style
 
 
 __all__ = ['Colormap3D', 'ColormapCubehelix', 'ColormapPerception', 'ColormapHLS', 'ColormapClassic',
-           'ColormapTransparent', 'cmaps', 'CMAP_CIRCULAR_DEFAULT', 'interpolate_color']
+           'ColormapTransparent', 'cmaps', 'interpolate_color']
 
 _log = logging.getLogger(__name__)
 
@@ -65,7 +64,7 @@ class Colormap3D(colors.Colormap, metaclass=abc.ABCMeta):
 
         """
         self._log.debug('Calling rgb_from_vector')
-        x, y, z = np.asarray(vector)
+        x, y, z = vector
         R = np.sqrt(x ** 2 + y ** 2 + z ** 2)
         R_max = vmax if vmax is not None else R.max() + 1E-30
         # FIRST color dimension: HUE (1D ring/angular direction)
@@ -498,14 +497,24 @@ def interpolate_color(fraction, start, end):
     return r1, r2, r3
 
 
-cmaps = {'cubehelix_standard': ColormapCubehelix(),
-         'cubehelix_reverse': ColormapCubehelix(reverse=True),
-         'cubehelix_circular': ColormapCubehelix(start=1, rot=1, minLight=0.5, maxLight=0.5, sat=2),
-         'perception_circular': ColormapPerception(),
-         'hls_circular': ColormapHLS(),
-         'classic_circular': ColormapClassic(),
-         'transparent_black': ColormapTransparent(0, 0, 0, [0, 1.]),
-         'transparent_white': ColormapTransparent(1, 1, 1, [0, 1.]),
-         'transparent_confidence': ColormapTransparent(0.2, 0.3, 0.2, [0.75, 0.])}
+class CMapNamespace(object):
 
-CMAP_CIRCULAR_DEFAULT = cmaps['cubehelix_circular']
+    def __init__(self):
+        self.cubehelix = ColormapCubehelix()
+        self.cubehelix_r = ColormapCubehelix(reverse=True)
+        self.cyclic_cubehelix = ColormapCubehelix(start=1, rot=1, minLight=0.5, maxLight=0.5, sat=2)
+        self.cyclic_perception = ColormapPerception()
+        self.cyclic_hls = ColormapHLS()
+        self.cyclic_classic = ColormapClassic()
+        self.transparent_black = ColormapTransparent(0, 0, 0, [0, 1.])
+        self.transparent_white = ColormapTransparent(1, 1, 1, [0, 1.])
+        self.transparent_confidence = ColormapTransparent(0.2, 0.3, 0.2, [0.75, 0.])
+
+    def __getitem__(self, key):
+        return self.__dict__[key]
+
+    def add_cmap_dict(self, **kwargs):
+        self.__dict__.update(kwargs)
+
+
+cmaps = CMapNamespace()

src/empyre/vis/decorators.py

@@ -94,8 +94,8 @@ def colorwheel(axis=None, cmap=None, ax_size='20%', loc='upper right', **kwargs)
     axis : :class:`~matplotlib.axes.Axes`, optional
         Axis to which the colorwheel is added, by default None, which will pick the last used axis via `gca`.
     cmap : str or `matplotlib.colors.Colormap`, optional
-        The Colormap that should be used for the colorwheel, defaults to `None`, which chooses a map specified by
-        `.colors.CMAP_CIRCULAR_DEFAULT`. Needs to be a :class:`~.colors.Colormap3D`.
+        The Colormap that should be used for the colorwheel, defaults to `None`, which chooses the
+        `.colors.cmaps.cyclic_cubehelix` colormap. Needs to be a :class:`~.colors.Colormap3D` to work correctly.
     ax_size : str or float, optional
         String or float determining the size of the inset axis used, defaults to `20%`.
     loc : str or pair of floats, optional
@@ -117,7 +117,7 @@ def colorwheel(axis=None, cmap=None, ax_size='20%', loc='upper right', **kwargs)
     ins_axes = inset_axes(axis, width=ax_size, height=ax_size, loc=loc)
     ins_axes.axis('off')
     if cmap is None:
-        cmap = colors.CMAP_CIRCULAR_DEFAULT
+        cmap = colors.cmaps.cyclic_cubehelix
     plt.sca(axis)  # Set focus back to parent axis!
     return cmap.plot_colorwheel(axis=ins_axes, **kwargs)
 
@@ -161,7 +161,7 @@ def quiverkey(quiv, field, axis=None, unit='', loc='lower right', **kwargs):
     quiv : Quiver instance
         The quiver instance returned by a call to quiver.
     field : `Field` or ndarray
-        The vector data as a `Field` or a numpy array (in the latter case, `vector=True` and `scale=1` are assumed).
+        The vector data as a `Field` or a numpy array (in the latter case, `vector=True` and `scale=1.0` are assumed).
     axis : :class:`~matplotlib.axes.AxesSubplot`, optional
         Axis to which the quiverkey is added, by default None, which will pick the last used axis via `gca`.
     unit: str, optional
@@ -183,7 +183,7 @@ def quiverkey(quiv, field, axis=None, unit='', loc='lower right', **kwargs):
     if axis is None:  # If no axis is set, find the current or create a new one:
         axis = plt.gca()
     if not isinstance(field, Field):  # Try to convert input to Field if it is not already one:
-        field = Field(data=np.asarray(field), scale=1, vector=True)
+        field = Field(data=np.asarray(field), scale=1.0, vector=True)
     length = field.amp.data.max()
     shift = 1 / field.squeeze().dim[1]  # equivalent to one pixel distance in axis coords!
     label = f'{length:.3g} {unit}'

src/empyre/vis/mplstyles/empyre-image.mplstyle

 ### MATPLOTLIB STYLESHEET FOR EMPYRE IMAGES
 
-text.usetex                   : True      ## use TeX to render text
-
 font.family                   : serif     ## default font family (use serifs)
 font.serif                    : cm        ## Computer Modern (LaTeX font)
 

src/empyre/vis/mplstyles/empyre-plot.mplstyle

 ### MATPLOTLIB STYLESHEET FOR EMPYRE PLOTS
 
-text.usetex                   : True      ## use TeX to render text
-
 font.family                   : serif     ## default font family (use serifs)
 font.serif                    : cm        ## Computer Modern (LaTeX font)