Merge commit '058ef52c' into demos-added

.gitignore

@@ -11,3 +11,5 @@ build/*
 docs/_build/*
 dist/*
 src/empyre/version.py
+htmlcov/
+coverage.xml

setup.cfg

@@ -7,7 +7,7 @@
 
 [metadata]
 name = empyre
-version = 0.2.0
+version = 0.2.1
 author = Jan Caron
 author-email = j.caron@fz-juelich.de
 description = Electron Microscopy Python Reconstruction
@@ -33,7 +33,7 @@ include_package_data = True
 package_dir =
     =src
 packages = find:
-python_requires = >=3.7
+python_requires = >=3.6
 setup_requires =
     setuptools
 install_requires =

src/empyre/fields/field.py

@@ -160,7 +160,7 @@ class Field(NDArrayOperatorsMixin):
 
     def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
         self._log.debug('Calling __array_ufunc__')
-        outputs = kwargs.pop('out', ())
+        # 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])
         # Cannot handle items that have __array_ufunc__ (other than our own).
@@ -178,22 +178,29 @@ class Field(NDArrayOperatorsMixin):
         # 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:
-            self._log.debug(f'__array_ufunc__ inputs: {len(inputs)}')
             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 = numeric.normalize_axis_tuple(axis, field.ndim)  # Takes care of pot. neg. indices, ensures tuple!
+                # axis=-1 reduces over the vector components, if they exist
+                # Takes care of pot. neg. indices, ensures tuple!
+                axis = numeric.normalize_axis_tuple(axis, len(field.shape))
                 kwargs['axis'] = axis  # Put normalized axis back into kwargs!
-                if axis in (ax_full, ax_full_wc):  # Full reduction (or reduction to just components) takes place:
-                    full_reduction = True
+                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!
@@ -210,7 +217,8 @@ class Field(NDArrayOperatorsMixin):
                 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}!'
-                assert all(scale == scale_new for scale in scales), err_msg
+                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:
@@ -218,7 +226,8 @@ class Field(NDArrayOperatorsMixin):
             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}!'
-                assert all(ncomp == ncomps[0] for ncomp in ncomps), err_msg
+                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)
@@ -408,34 +417,35 @@ class Field(NDArrayOperatorsMixin):
             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]
-        if axis == 'z':  # TODO: Somehow needed... don't know why, maybe because origin='lower' in imshow?
-            angle *= -1
+            # 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 = interpolation.rotate(self.data, angle, axes=axes, **kwargs)
+            data_new = interpolation.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([interpolation.rotate(c, angle, axes=axes, **kwargs) for c in comps], axis=-1)
+                data_new = np.stack([interpolation.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 = [i for i in (0, 1, 2) if i not in axes][0]  # [0] because list!
+                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):
-                # TODO: Again a minus sign needed before levi_civita, still not sure why (see angle above)!
-                quat_axis = tuple([-levi_civita if i == rot_axis else 0 for i in (2, 1, 0)])
+                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))
-                data_new = quat.matrix.dot(data_new.reshape((-1, 3)).T).T.reshape(self.shape)  # T needed b.c. ordering!
+                # 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 = interpolation.rotate(u_comp, angle, axes=axes, **kwargs)
-                v_rot = interpolation.rotate(v_comp, angle, axes=axes, **kwargs)
+                u_rot = interpolation.rotate(u_comp, np_angle, axes=axes, **kwargs)
+                v_rot = interpolation.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:
-                # TODO: Again a minus sign needed for vector rotation, still not sure why (see angle above)!
-                ang_rad = np.deg2rad(-angle)
+                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)
@@ -470,29 +480,33 @@ class Field(NDArrayOperatorsMixin):
         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)
+            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)
+                    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_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_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:
@@ -662,6 +676,27 @@ class Field(NDArrayOperatorsMixin):
         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 = numeric.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.
 
@@ -764,21 +799,26 @@ class Field(NDArrayOperatorsMixin):
 
         """
         self._log.debug('Calling clip')
-        amp = self.amp.data
+        # 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 = (amp - amp.mean()) < (sigma * amp.std())
+            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:
-        amp_masked = np.where(mask, amp, np.nan)
+        indicator_masked = np.where(mask, indicator, np.nan)
         if vmin is None:
-            vmin = np.nanmin(amp_masked)
+            vmin = np.nanmin(indicator_masked)
         if vmax is None:
-            vmax = np.nanmax(amp_masked)
-        if self.vector:  # Vector fields need to scale components according to masked amplitude
-            mask_vec = (amp <= vmax)[..., None]  # Only vmax is important for vectors! mask_vec broadcast to components!
-            data = np.where(mask_vec, self.data, vmax * self.data/amp[..., None])  # Scale outliers to vmax!
+            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 = np.clip(self.data, vmin, vmax)
-        return Field(data, self.scale, self.vector)
+            data_new = np.clip(self.data, vmin, vmax)
+        return Field(data_new, self.scale, self.vector)

src/empyre/io/field_plugins/vtk.py

@@ -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,10 @@ 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!'
         data = interp_to_regular_grid(point_array, data_array, scale, **kwargs)
     else:
         raise TypeError('Data type of {} not understood!'.format(output))

src/empyre/utils/misc.py

@@ -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
     -------
@@ -67,12 +73,12 @@ def interp_to_regular_grid(points, values, scale, scale_factor=1, step=1, convex
     _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)))
-    x = x_min + scale[2] * (np.arange(dim[2]) + 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!
+    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:
-    test = np.asarray(list(itertools.product(z, y, x)))
-    points_euc = np.fliplr(test)
+    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:
@@ -93,12 +99,25 @@ 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)

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)

tests/conftest.py

+import os
+
+import pytest
+import numpy as np
+
+from empyre.fields import Field
+
+
+@pytest.fixture
+def fielddata_path():
+    return os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_fielddata')
+
+
+@pytest.fixture
+def vector_data():
+    magnitude = np.zeros((4, 4, 4, 3))
+    magnitude[1:-1, 1:-1, 1:-1] = 1
+    return Field(magnitude, 10.0, vector=True)
+
+
+@pytest.fixture
+def vector_data_asymm():
+    shape = (5, 7, 11, 3)
+    data = np.linspace(0, 1, np.prod(shape))
+    return Field(data.reshape(shape), 10.0, vector=True)
+
+
+@pytest.fixture
+def vector_data_asymm_2d():
+    shape = (5, 7, 2)
+    data = np.linspace(0, 1, np.prod(shape))
+    return Field(data.reshape(shape), 10.0, vector=True)
+
+
+@pytest.fixture
+def vector_data_asymmcube():
+    shape = (3, 3, 3, 3)
+    data = np.linspace(0, 1, np.prod(shape))
+    return Field(data.reshape(shape), 10.0, vector=True)
+
+
+@pytest.fixture
+def scalar_data():
+    magnitude = np.zeros((4, 4, 4))
+    magnitude[1:-1, 1:-1, 1:-1] = 1
+    return Field(magnitude, 10.0, vector=False)
+
+
+@pytest.fixture
+def scalar_data_asymm():
+    shape = (5, 7, 2)
+    data = np.linspace(0, 1, np.prod(shape))
+    return Field(data.reshape(shape), 10.0, vector=False)

tests/fields/test_fielddata.py

+# -*- coding: utf-8 -*-
+"""Testcase for the magdata module."""
+
+import pytest
+
+from numbers import Number
+
+import numpy as np
+import numpy.testing
+
+from empyre.fields import Field
+
+from utils import assert_allclose
+
+
+def test_copy(vector_data):
+    vector_data = vector_data.copy()
+    # Make sure it is a new object
+    assert vector_data != vector_data, 'Unexpected behaviour in copy()!'
+    assert np.allclose(vector_data, vector_data)
+
+
+def test_bin(vector_data):
+    binned_data = vector_data.bin(2)
+    reference = 1 / 8. * np.ones((2, 2, 2, 3))
+    assert_allclose(binned_data, reference,
+                    err_msg='Unexpected behavior in scale_down()!')
+    assert_allclose(binned_data.scale, (20, 20, 20),
+                    err_msg='Unexpected behavior in scale_down()!')
+
+
+def test_zoom(vector_data):
+    zoomed_test = vector_data.zoom(2, order=0)
+    reference = np.zeros((8, 8, 8, 3))
+    reference[2:6, 2:6, 2:6] = 1
+    assert_allclose(zoomed_test, reference,
+                    err_msg='Unexpected behavior in zoom()!')
+    assert_allclose(zoomed_test.scale, (5, 5, 5),
+                    err_msg='Unexpected behavior in zoom()!')
+
+
+@pytest.mark.parametrize(
+    'mode', [
+        'constant',
+        'edge',
+        'wrap'
+    ]
+)
+@pytest.mark.parametrize(
+    'pad_width,np_pad', [
+        (1, ((1, 1), (1, 1), (1, 1), (0, 0))),
+        ((1, 2, 3), ((1, 1), (2, 2), (3, 3), (0, 0))),
+        (((1, 2), (3, 4), (5, 6)), ((1, 2), (3, 4), (5, 6), (0, 0)))
+    ]
+)
+def test_pad(vector_data, mode, pad_width, np_pad):
+    magdata_test = vector_data.pad(pad_width, mode=mode)
+    reference = np.pad(vector_data, np_pad, mode=mode)
+    assert_allclose(magdata_test, reference,
+                    err_msg='Unexpected behavior in pad()!')
+
+
+@pytest.mark.parametrize(
+    'axis', [-1, 3]
+)
+def test_component_reduction(vector_data, axis):
+    # axis=-1 is supposed to reduce over the component dimension, if it exists. axis=3 should do the same here!
+    res = np.sum(vector_data, axis=axis)
+    ref = np.zeros((4, 4, 4))
+    ref[1:-1, 1:-1, 1:-1] = 3
+    assert res.shape == ref.shape, 'Shape mismatch!'
+    assert_allclose(res, ref, err_msg="Unexpected behavior of axis keyword")
+    assert isinstance(res, Field), 'Result is not a Field object!'
+    assert not res.vector, 'Result is a vector field, but should be reduced to a scalar!'
+
+
+@pytest.mark.parametrize(
+    'axis', [(0, 1, 2), (2, 1, 0), None, (-4, -3, -2)]
+)
+def test_full_reduction(vector_data, axis):
+    res = np.sum(vector_data, axis=axis)
+    ref = np.zeros((3,))
+    ref[:] = 8
+    assert res.shape == ref.shape
+    assert_allclose(res, ref, err_msg="Unexpected behavior of full or default reduction")
+    assert isinstance(res, np.ndarray)
+
+
+@pytest.mark.parametrize(
+    'axis', [-1, 2]
+)
+def test_last_reduction_scalar(scalar_data, axis):
+    # axis=-1 is supposed to reduce over the component dimension if it exists.
+    # In this case it doesn't!
+    res = np.sum(scalar_data, axis=axis)
+    ref = np.zeros((4, 4))
+    ref[1:-1, 1:-1] = 2
+    assert res.shape == ref.shape
+    assert_allclose(res, ref, err_msg="Unexpected behavior of axis keyword")
+    assert isinstance(res, Field)
+    assert not res.vector
+
+
+@pytest.mark.parametrize(
+    'axis', [(0, 1, 2), (2, 1, 0), None, (-1, -2, -3)]
+)
+def test_full_reduction_scalar(scalar_data, axis):
+    res = np.sum(scalar_data, axis=axis)
+    ref = 8
+    assert res.shape == ()
+    assert_allclose(res, ref, err_msg="Unexpected behavior of full or default reduction")
+    assert isinstance(res, Number)
+
+
+def test_binary_operator_vector_number(vector_data):
+    res = vector_data + 1
+    ref = np.ones((4, 4, 4, 3))
+    ref[1:-1, 1:-1, 1:-1] = 2
+    assert res.shape == ref.shape
+    assert_allclose(res, ref, err_msg="Unexpected behavior of addition")
+    assert isinstance(res, Field)
+    assert res.vector
+
+
+def test_binary_operator_vector_scalar(vector_data, scalar_data):
+    res = vector_data + scalar_data
+    ref = np.zeros((4, 4, 4, 3))
+    ref[1:-1, 1:-1, 1:-1] = 2
+    assert res.shape == ref.shape
+    assert_allclose(res, ref, err_msg="Unexpected behavior of addition")
+    assert isinstance(res, Field)
+    assert res.vector
+
+
+def test_binary_operator_vector_vector(vector_data):
+    res = vector_data + vector_data
+    ref = np.zeros((4, 4, 4, 3))
+    ref[1:-1, 1:-1, 1:-1] = 2
+    assert res.shape == ref.shape
+    assert_allclose(res, ref, err_msg="Unexpected behavior of addition")
+    assert isinstance(res, Field)
+    assert res.vector
+
+
+@pytest.mark.xfail
+def test_binary_operator_vector_broadcast(vector_data):
+    # Broadcasting between vector fields is currently not implemented
+    second = np.zeros((4, 4, 3))
+    second[1:-1, 1:-1] = 1
+    second = Field(second, 10.0, vector=True)
+    res = vector_data + second
+    ref = np.zeros((4, 4, 4, 3))
+    ref[1:-1, 1:-1, 1:-1] = 1
+    ref[:, 1:-1, 1:-1] += 1
+    assert res.shape == ref.shape
+    assert_allclose(res, ref, err_msg="Unexpected behavior of addition")
+    assert isinstance(res, Field)
+    assert res.vector
+
+
+def test_mask(vector_data):
+    mask = vector_data.mask
+    reference = np.zeros((4, 4, 4))
+    reference[1:-1, 1:-1, 1:-1] = True
+    assert_allclose(mask, reference,
+                    err_msg='Unexpected behavior in mask attribute!')
+
+
+def test_get_vector(vector_data):
+    mask = vector_data.mask
+    vector = vector_data.get_vector(mask)
+    reference = np.ones(np.sum(mask) * 3)
+    assert_allclose(vector, reference,
+                    err_msg='Unexpected behavior in get_vector()!')
+
+
+def test_set_vector(vector_data):
+    mask = vector_data.mask
+    vector = 2 * np.ones(np.sum(mask) * 3)
+    vector_data.set_vector(vector, mask)
+    reference = np.zeros((4, 4, 4, 3))
+    reference[1:-1, 1:-1, 1:-1] = 2
+    assert_allclose(vector_data, reference,
+                    err_msg='Unexpected behavior in set_vector()!')
+
+
+def test_flip(vector_data_asymm):
+    field_flipz = vector_data_asymm.flip(0)
+    field_flipy = vector_data_asymm.flip(1)
+    field_flipx = vector_data_asymm.flip(2)
+    field_flipxy = vector_data_asymm.flip((1, 2))
+    field_flipdefault = vector_data_asymm.flip()
+    field_flipcomp = vector_data_asymm.flip(-1)
+    assert_allclose(np.flip(vector_data_asymm.data, axis=0) * [1, 1, -1], field_flipz.data,
+                    err_msg='Unexpected behavior in flip()! (z)')
+    assert_allclose(np.flip(vector_data_asymm.data, axis=1) * [1, -1, 1], field_flipy.data,
+                    err_msg='Unexpected behavior in flip()! (y)')
+    assert_allclose(np.flip(vector_data_asymm.data, axis=2) * [-1, 1, 1], field_flipx.data,
+                    err_msg='Unexpected behavior in flip()! (x)')
+    assert_allclose(np.flip(vector_data_asymm.data, axis=(1, 2)) * [-1, -1, 1], field_flipxy.data,
+                    err_msg='Unexpected behavior in flip()! (xy)')
+    assert_allclose(np.flip(vector_data_asymm.data, axis=(0, 1, 2)) * [-1, -1, -1], field_flipdefault.data,
+                    err_msg='Unexpected behavior in flip()! (default)')
+    assert_allclose(np.flip(vector_data_asymm.data, axis=-1) * [1, 1, 1], field_flipcomp.data,
+                    err_msg='Unexpected behavior in flip()! (components)')
+
+
+def test_unknown_num_of_components():
+    shape = (5, 7, 7)
+    data = np.linspace(0, 1, np.prod(shape))
+    with pytest.raises(AssertionError):
+        Field(data.reshape(shape), 10.0, vector=True)
+
+
+def test_repr(vector_data_asymm):
+    string_repr = repr(vector_data_asymm)
+    data_str = str(vector_data_asymm.data)
+    string_ref = f'Field(data={data_str}, scale=(10.0, 10.0, 10.0), vector=True)'
+    print(f'reference: {string_ref}')
+    print(f'repr output: {string_repr}')
+    assert string_repr == string_ref, 'Unexpected behavior in __repr__()!'
+
+
+def test_str(vector_data_asymm):
+    string_str = str(vector_data_asymm)
+    string_ref = 'Field(dim=(5, 7, 11), scale=(10.0, 10.0, 10.0), vector=True, ncomp=3)'
+    print(f'reference: {string_str}')
+    print(f'str output: {string_str}')
+    assert string_str == string_ref, 'Unexpected behavior in __str__()!'
+
+
+@pytest.mark.parametrize(
+    "index,t,scale", [
+        ((0, 1, 2), tuple, None),
+        ((0, ), Field, (2., 3.)),
+        (0, Field, (2., 3.)),
+        ((0, 1, 2, 0), float, None),
+        ((0, 1, 2, 0), float, None),
+        ((..., 0), Field, (1., 2., 3.)),
+        ((0, slice(1, 3), 2), Field, (2.,)),
+    ]
+)
+def test_getitem(vector_data, index, t, scale):
+    vector_data.scale = (1., 2., 3.)
+    data_index = index
+    res = vector_data[index]
+    assert_allclose(res, vector_data.data[data_index])
+    assert isinstance(res, t)
+    if t is Field:
+        assert res.scale == scale
+
+
+def test_from_scalar_field(scalar_data):
+    sca_x, sca_y, sca_z = [i * scalar_data for i in range(1, 4)]
+    field_comb = Field.from_scalar_fields([sca_x, sca_y, sca_z])
+    assert field_comb.vector
+    assert field_comb.scale == scalar_data.scale
+    assert_allclose(sca_x, field_comb.comp[0])
+    assert_allclose(sca_y, field_comb.comp[1])
+    assert_allclose(sca_z, field_comb.comp[2])
+
+
+def test_squeeze():
+    magnitude = np.zeros((4, 1, 4, 3))
+    field = Field(magnitude, (1., 2., 3.), vector=True)
+    sq = field.squeeze()
+    assert sq.shape == (4, 4, 3)
+    assert sq.dim == (4, 4)
+    assert sq.scale == (1., 3.)
+
+
+def test_gradient():
+    pass
+
+
+def test_gradient_1d():
+    pass
+
+
+def test_curl():
+    pass
+
+
+def test_curl_2d():
+    pass
+
+
+def test_clip_scalar_noop():
+    shape = (3, 3, 3)
+    data = np.linspace(-2, 1, np.prod(shape)).reshape(shape)
+    field = Field(data, (1., 2., 3.), vector=False)
+    assert_allclose(field, field.clip())
+
+
+def test_clip_scalar_minmax():
+    shape = (3, 3, 3)
+    data = np.linspace(-2, 1, np.prod(shape)).reshape(shape)
+    field = Field(data, (1., 2., 3.), vector=False)
+    assert_allclose(np.clip(data, -1, 0.1), field.clip(vmin=-1, vmax=0.1))
+
+
+def test_clip_scalar_sigma():
+    shape = (3, 3, 3)
+    data = np.linspace(-2, 1, np.prod(shape)).reshape(shape)
+    data[0, 0, 0] = 1e6
+    field = Field(data, (1., 2., 3.), vector=False)
+    # We clip off the one outlier
+    assert_allclose(np.clip(data, -2, 1), field.clip(sigma=5))
+    assert field.clip(sigma=5)[0, 0, 0] == 1
+
+
+def test_clip_scalar_mask():
+    shape = (3, 3, 3)
+    data = np.linspace(-2, 1, np.prod(shape)).reshape(shape)
+    mask = np.zeros(shape, dtype=bool)
+    mask[0, 0, 0] = True
+    mask[0, 0, 1] = True
+    field = Field(data, (1., 2., 3.), vector=False)
+    assert_allclose(np.clip(data, data[0, 0, 0], data[0, 0, 1]), field.clip(mask=mask))
+
+
+def test_clip_vector_noop():
+    shape = (3, 3, 3, 3)
+    data = np.linspace(-2, 1, np.prod(shape)).reshape(shape)
+    field = Field(data, (1., 2., 3.), vector=True)
+    assert_allclose(field, field.clip())
+
+
+def test_clip_vector_max():
+    shape = (3, 3, 3, 3)
+    data = np.linspace(-2, 1, np.prod(shape)).reshape(shape)
+    field = Field(data, (1., 2., 3.), vector=True)
+    res = field.clip(vmax=0.1)
+    assert_allclose(np.max(res.amp), 0.1)
+
+
+def test_clip_vector_sigma():
+    shape = (3, 3, 3, 3)
+    data = np.linspace(-2, 1, np.prod(shape)).reshape(shape)
+    data[0, 0, 0] = (1e6, 1e6, 1e6)
+    field = Field(data, (1., 2., 3.), vector=True)
+    # We clip off the one outlier
+    res = field.clip(sigma=5)
+    assert np.max(res.amp) < 1e3
+
+
+# TODO: HyperSpy would need to be installed for the following tests (slow...):
+# def test_from_signal()
+#     raise NotImplementedError()
+#
+# def test_to_signal()
+#     raise NotImplementedError()

tests/fields/test_rotate.py

+import pytest
+
+from utils import assert_allclose
+
+from empyre.fields import Field
+import numpy as np
+
+
+@pytest.mark.parametrize(
+    'axis', ['x', 'y', 'z']
+)
+def test_rot90_360(vector_data_asymm, axis):
+    assert_allclose(vector_data_asymm.rot90(axis=axis).rot90(axis=axis).rot90(axis=axis).rot90(axis=axis),
+                    vector_data_asymm,
+                    err_msg=f'Unexpected behavior in rot90()! {axis}')
+
+
+@pytest.mark.parametrize(
+    'rot_axis,flip_axes', [
+        ('x', (0, 1)),
+        ('y', (0, 2)),
+        ('z', (1, 2))
+    ]
+)
+def test_rot90_180(vector_data_asymm, rot_axis, flip_axes):
+    res = vector_data_asymm.rot90(axis=rot_axis).rot90(axis=rot_axis)
+    ref = vector_data_asymm.flip(axis=flip_axes)
+    assert_allclose(res, ref, err_msg=f'Unexpected behavior in rot90()! {rot_axis}')
+
+
+@pytest.mark.parametrize(
+    'rot_axis', [
+        'x',
+        'y',
+        'z',
+    ]
+)
+def test_rotate_compare_rot90_1(vector_data_asymmcube, rot_axis):
+    res = vector_data_asymmcube.rotate(angle=90, axis=rot_axis)
+    ref = vector_data_asymmcube.rot90(axis=rot_axis)
+    print("input", vector_data_asymmcube.data)
+    print("ref", res.data)
+    print("res", ref.data)
+    assert_allclose(res, ref, err_msg=f'Unexpected behavior in rotate()! {rot_axis}')
+
+
+def test_rot90_manual():
+    data = np.zeros((3, 3, 3, 3))
+    diag = np.array((1, 1, 1))
+    diag_unity = diag / np.sqrt(np.sum(diag**2))
+    data[0, 0, 0] = diag_unity
+    data = Field(data, 10, vector=True)
+    print("data", data.data)
+
+    rot90_x = np.zeros((3, 3, 3, 3))
+    # Axis order z, y, x; vector components x, y, z
+    rot90_x[0, 2, 0] = diag_unity * (1, -1, 1)
+    rot90_x = Field(rot90_x, 10, vector=True)
+    print("rot90_x", rot90_x.data)
+    print("data rot90 x", data.rot90(axis='x').data)
+
+    rot90_y = np.zeros((3, 3, 3, 3))
+    # Axis order z, y, x; vector components x, y, z
+    rot90_y[2, 0, 0] = diag_unity * (1, 1, -1)
+    rot90_y = Field(rot90_y, 10, vector=True)
+    print("rot90_y", rot90_y.data)
+    print("data rot90 y", data.rot90(axis='y').data)
+
+    rot90_z = np.zeros((3, 3, 3, 3))
+    # Axis order z, y, x; vector components x, y, z
+    rot90_z[0, 0, 2] = diag_unity * (-1, 1, 1)
+    rot90_z = Field(rot90_z, 10, vector=True)
+    print("rot90_z", rot90_z.data)
+    print("data rot90 z", data.rot90(axis='z').data)
+
+    assert_allclose(rot90_x, data.rot90(axis='x'), err_msg='Unexpected behavior in rot90("x")!')
+    assert_allclose(rot90_y, data.rot90(axis='y'), err_msg='Unexpected behavior in rot90("y")!')
+    assert_allclose(rot90_z, data.rot90(axis='z'), err_msg='Unexpected behavior in rot90("z")!')
+
+
+def test_rot45_manual():
+    data = np.zeros((3, 3, 3, 3))
+    data[0, 0, 0] = (1, 1, 1)
+    data = Field(data, 10, vector=True)
+    print("data", data.data)
+
+    rot45_x = np.zeros((3, 3, 3, 3))
+    # Axis order z, y, x; vector components x, y, z
+    rot45_x[0, 1, 0] = (1, 0, np.sqrt(2))
+    rot45_x = Field(rot45_x, 10, vector=True)
+    print("rot45_x", rot45_x.data)
+    # Disable spline interpolation, use nearest instead
+    res_rot45_x = data.rotate(45, axis='x', order=0)
+    print("data rot45 x", res_rot45_x.data)
+
+    rot45_y = np.zeros((3, 3, 3, 3))
+    # Axis order z, y, x; vector components x, y, z
+    rot45_y[1, 0, 0] = (np.sqrt(2), 1, 0)
+    rot45_y = Field(rot45_y, 10, vector=True)
+    print("rot45_y", rot45_y.data)
+    # Disable spline interpolation, use nearest instead
+    res_rot45_y = data.rotate(45, axis='y', order=0)
+    print("data rot45 y", res_rot45_y.data)
+
+    rot45_z = np.zeros((3, 3, 3, 3))
+    # Axis order z, y, x; vector components x, y, z
+    rot45_z[0, 0, 1] = (0, np.sqrt(2), 1)
+    rot45_z = Field(rot45_z, 10, vector=True)
+    print("rot45_z", rot45_z.data)
+    # Disable spline interpolation, use nearest instead
+    res_rot45_z = data.rotate(45, axis='z', order=0)
+    print("data rot45 z", res_rot45_z.data)
+
+    assert_allclose(rot45_x, res_rot45_x, err_msg='Unexpected behavior in rotate(45, "x")!')
+    assert_allclose(rot45_y, res_rot45_y, err_msg='Unexpected behavior in rotate(45, "y")!')
+    assert_allclose(rot45_z, res_rot45_z, err_msg='Unexpected behavior in rotate(45, "z")!')
+
+
+def test_rot90_2d_360(vector_data_asymm_2d):
+    assert_allclose(vector_data_asymm_2d.rot90().rot90().rot90().rot90(), vector_data_asymm_2d,
+                    err_msg='Unexpected behavior in 2D rot90()!')
+
+
+def test_rot90_2d_180(vector_data_asymm_2d):
+    res = vector_data_asymm_2d.rot90().rot90()
+    ref = vector_data_asymm_2d.flip()
+    assert_allclose(res, ref, err_msg='Unexpected behavior in 2D rot90()!')
+
+
+@pytest.mark.parametrize(
+    'k', [0, 1, 2, 3, 4]
+)
+def test_rot90_comp_2d_with_3d(vector_data_asymm_2d, k):
+    data_x, data_y = [comp.data[np.newaxis, :, :] for comp in vector_data_asymm_2d.comp]
+    data_z = np.zeros_like(data_x)
+    data_3d = np.stack([data_x, data_y, data_z], axis=-1)
+    vector_data_asymm_3d = Field(data_3d, scale=10, vector=True)
+    print(f'2D shape, scale: {vector_data_asymm_2d.shape, vector_data_asymm_2d.scale}')
+    print(f'3D shape, scale: {vector_data_asymm_3d.shape, vector_data_asymm_3d.scale}')
+    vector_data_rot_2d = vector_data_asymm_2d.rot90(k=k)
+    vector_data_rot_3d = vector_data_asymm_3d.rot90(k=k, axis='z')
+    print(f'2D shape after rot: {vector_data_rot_2d.shape}')
+    print(f'3D shape after rot: {vector_data_rot_3d.shape}')
+    assert_allclose(vector_data_rot_2d, vector_data_rot_3d[0, :, :, :2], err_msg='Unexpected behavior in 2D rot90()!')
+
+
+@pytest.mark.parametrize(
+    'angle', [90, 45, 23, 11.5]
+)
+def test_rotate_comp_2d_with_3d(vector_data_asymm_2d, angle):
+    data_x, data_y = [comp.data[np.newaxis, :, :] for comp in vector_data_asymm_2d.comp]
+    data_z = np.zeros_like(data_x)
+    data_3d = np.stack([data_x, data_y, data_z], axis=-1)
+    vector_data_asymm_3d = Field(data_3d, scale=10, vector=True)
+    print(f'2D shape, scale: {vector_data_asymm_2d.shape, vector_data_asymm_2d.scale}')
+    print(f'3D shape, scale: {vector_data_asymm_3d.shape, vector_data_asymm_3d.scale}')
+    r2d = vector_data_asymm_2d.rotate(angle)
+    r3d = vector_data_asymm_3d.rotate(angle, axis='z')
+    print(f'2D shape after rot: {r2d.shape}')
+    print(f'3D shape after rot: {r3d.shape}')
+    assert_allclose(r2d, r3d[0, :, :, :2], err_msg='Unexpected behavior in 2D rotate()!')
+
+
+@pytest.mark.parametrize(
+    'angle', [180, 360, 90, 45, 23, 11.5],
+)
+@pytest.mark.parametrize(
+    'axis', ['x', 'y', 'z'],
+)
+def test_rotate_scalar(vector_data_asymm, angle, axis):
+
+    data = np.zeros((1, 2, 2, 3))
+    data[0, 0, 0] = 1
+    field = Field(data, scale=10., vector=True)
+    print(field)
+    print(field.amp)
+
+    assert_allclose(
+        field.rotate(angle, axis=axis).amp,
+        field.amp.rotate(angle, axis=axis)
+    )
+
+
+@pytest.mark.parametrize(
+    'angle,order', [(180, 3), (360, 3), (90, 3), (45, 0), (23, 0), (11.5, 0)],
+)
+@pytest.mark.parametrize(
+    'axis', ['x', 'y', 'z'],
+)
+@pytest.mark.parametrize(
+    'reshape', [True, False],
+)
+def test_rotate_scalar_asymm(vector_data_asymm, angle, axis, order, reshape):
+    assert_allclose(
+        vector_data_asymm.rotate(angle, axis=axis, reshape=reshape, order=order).amp,
+        vector_data_asymm.amp.rotate(angle, axis=axis, reshape=reshape, order=order)
+    )
+
+
+@pytest.mark.parametrize(
+    'axis', ['x', 'y', 'z'],
+)
+@pytest.mark.parametrize(
+    'k', [0, 1, 2, 3, 4],
+)
+def test_rot90_scalar(vector_data_asymm, axis, k):
+    assert_allclose(
+        vector_data_asymm.amp.rot90(k=k, axis=axis),
+        vector_data_asymm.rot90(k=k, axis=axis).amp
+    )

tests/test_fields/test_field.py

-# -*- coding: utf-8 -*-
-"""Testcase for the magdata module."""
-
-# import os
-import unittest
-
-# import numpy as np
-# from numpy.testing import assert_allclose
-
-# from empyre.fields.field import Field
-
-
-class TestCaseField(unittest.TestCase):
-
-    def test_stupid(self):
-        print('FOUND')
-        assert True
-
-
-def test_simple():
-    print('METHOD')
-    assert True

tests/utils.py

+import numpy
+import numpy.testing
+
+
+def assert_allclose(actual, desired, rtol=1e-07, atol=1e-08, equal_nan=True, err_msg='', verbose=True):
+    return numpy.testing.assert_allclose(
+        actual=actual,
+        desired=desired,
+        rtol=rtol,
+        atol=atol,
+        equal_nan=equal_nan,
+        err_msg=err_msg,
+        verbose=verbose,
+    )