Merge pull request #453 from pc494/enh-strain-mapping-class

Strain Mapping Enhancement: StrainMap class, with basis changing method

pyxem/signals/__init__.py

@@ -48,6 +48,32 @@ def push_metadata_through(dummy, *args, **kwargs):
 
     return dummy, args, kwargs
 
+def transfer_signal_axes(new_signal, old_signal):
+    """ Transfers signal axis calibrations from an old signal to a new
+    signal produced from it by a method or a generator.
+
+    Parameters
+    ----------
+    new_signal : Signal
+        The product signal with undefined signal axes.
+    old_signal : Signal
+        The parent signal with calibrated signal axes.
+
+    Returns
+    -------
+    new_signal : Signal
+        The new signal with calibrated signal axes.
+    """
+
+    for i in range(old_signal.axes_manager.signal_dimension):
+        ax_new = new_signal.axes_manager.signal_axes[i]
+        ax_old = old_signal.axes_manager.signal_axes[i]
+        ax_new.name = ax_old.name
+        ax_new.scale = ax_old.scale
+        ax_new.units = ax_old.units
+
+    return new_signal
+
 
 def transfer_navigation_axes(new_signal, old_signal):
     """ Transfers navigation axis calibrations from an old signal to a new

pyxem/signals/strain_map.py

@@ -0,0 +1,124 @@
+# -*- coding: utf-8 -*-
+# Copyright 2017-2019 The pyXem developers
+#
+# This file is part of pyXem.
+#
+# pyXem is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# pyXem is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with pyXem.  If not, see <http://www.gnu.org/licenses/>.
+
+from hyperspy.signals import Signal2D
+import numpy as np
+from pyxem.signals import push_metadata_through, transfer_signal_axes
+
+
+def _get_rotation_matrix(x_new):
+    """Calculate the rotation matrix mapping [1,0] to x_new.
+
+    Parameters
+    ----------
+    x_new : list
+        Coordinates of a point on the new 'x' axis.
+     
+    Returns
+    -------
+    R : 2 x 2 numpy.array()
+        The rotation matrix.
+    """
+    try:
+        rotation_angle = np.arctan(x_new[1] / x_new[0])
+    except ZeroDivisionError:  # Taking x --> y
+        rotation_angle = np.deg2rad(90)
+
+    # angle sign agrees with https://en.wikipedia.org/wiki/Rotation_matrix
+    R = np.array([[np.cos(rotation_angle), -np.sin(rotation_angle)],
+                  [np.sin(rotation_angle), np.cos(rotation_angle)]])
+    return R
+
+
+class StrainMap(Signal2D):
+    """
+    Class for storing strain maps, if created within pyxem conventions are:
+    The 'y-axis' is 90 degrees from the 'x-axis'
+    Positive rotations are anticlockwise.
+    """
+
+    _signal_type = "strain_map"
+
+    def __init__(self, *args, **kwargs):
+        self, args, kwargs = push_metadata_through(self, *args, **kwargs)
+        super().__init__(*args, **kwargs)
+
+        # check init dimension are correct
+
+        if 'current_basis_x' in kwargs.keys():
+            self.current_basis_x = kwargs['current_basis_x']
+        else:
+            self.current_basis_x = [1, 0]
+
+        self.current_basis_y = np.matmul(np.asarray([[0, 1], [-1, 0]]), self.current_basis_x)
+
+    def rotate_strain_basis(self, x_new):
+        """ Rotates a strain map to a new basis.
+
+        Parameters
+        ----------
+        x_new : list
+            The coordinates of a point on the new 'x' axis
+
+        Returns
+        -------
+        StrainMap :
+            StrainMap in the new (rotated) basis.
+
+        Notes
+        -----
+        Conventions are described in the class documentation.
+        
+        We follow mathmatical formalism described in:
+        "https://www.continuummechanics.org/stressxforms.html" (August 2019)
+        """
+
+        def apply_rotation(transposed_strain_map, R):
+            """ Rotates a strain matrix to a new basis, for which R maps x_old to x_new """
+            sigmaxx_old = transposed_strain_map[0]
+            sigmayy_old = transposed_strain_map[1]
+            sigmaxy_old = transposed_strain_map[2]
+
+            z = np.asarray([[sigmaxx_old, sigmaxy_old],
+                            [sigmaxy_old, sigmayy_old]])
+
+            new = np.matmul(R.T, np.matmul(z, R))
+            return [new[0, 0], new[1, 1], new[0, 1], transposed_strain_map[3]]
+
+        def apply_rotation_complete(self, R):
+            """ Mapping solution to return a (unclassed) strain map in a new basis """
+            from hyperspy.api import transpose
+            transposed = transpose(self)[0]
+            transposed_to_new_basis = transposed.map(apply_rotation, R=R, inplace=False)
+            return transposed_to_new_basis.T
+
+        """ Core functionality """
+
+        if self.current_basis_x != [1, 0]:
+            # this takes us back to [1,0] if our current map is in a diferent basis
+            R = _get_rotation_matrix(self.current_basis_x).T
+            strain_map_core = apply_rotation_complete(self, R)
+        else:
+            strain_map_core = self
+
+        R = _get_rotation_matrix(x_new)
+        transposed_to_new_basis = apply_rotation_complete(strain_map_core, R)
+        meta_dict = self.metadata.as_dictionary()
+
+        strainmap = StrainMap(transposed_to_new_basis, current_basis_x=x_new, metadata=meta_dict)
+        return transfer_signal_axes(strainmap,self)

pyxem/signals/tensor_field.py

@@ -21,6 +21,7 @@
 from scipy.linalg import polar
 from hyperspy.utils import stack
 import math
+from pyxem.signals.strain_map import StrainMap
 
 """
 Signal class for Tensor Fields
@@ -111,4 +112,4 @@ def get_strain_maps(self):
 
         strain_results = stack([e11, e22, e12, theta])
 
-        return strain_results
+        return StrainMap(strain_results)

pyxem/tests/test_generators/test_displacement_gradient_tensor_generator.py

@@ -1,9 +1,4 @@
-from pyxem.generators.displacement_gradient_tensor_generator import \
-    get_DisplacementGradientMap, get_single_DisplacementGradientTensor
-import hyperspy.api as hs
-import pytest
-import numpy as np
-decimal = 2  # -*- coding: utf-8 -*-
+# -*- coding: utf-8 -*-
 # Copyright 2017-2019 The pyXem developers
 #
 # This file is part of pyXem.
@@ -21,6 +16,12 @@
 # You should have received a copy of the GNU General Public License
 # along with pyXem.  If not, see <http://www.gnu.org/licenses/>.
 
+from pyxem.generators.displacement_gradient_tensor_generator import \
+    get_DisplacementGradientMap, get_single_DisplacementGradientTensor
+import hyperspy.api as hs
+import pytest
+import numpy as np
+
 
 def vector_operation(z, M):
     """
@@ -142,37 +143,3 @@ def test_weight_function_behaviour():
     deformed = hs.signals.Signal2D(np.asarray([[vectors, vectors], [vectors, vectors]]))
     strain_map = get_DisplacementGradientMap(deformed, multi_vector_array, weights=weights).get_strain_maps()
     np.testing.assert_almost_equal(strain_map.inav[0].isig[0, 0].data[0], -1.0166666 + 1, decimal=2)
-
-
-""" These test will be operational once a basis change functionality is introduced """
-
-
-@pytest.mark.skip(reason="basis change functionality not yet implemented")
-def test_rotation(xy_vectors, right_handed, left_handed, multi_vector):  # pragma: no cover
-    """
-    We should always measure the same rotations, regardless of basis (as long as it's right handed)
-    """
-    xy_rot = xy_vectors.inav[3].data
-    rh_rot = right_handed.inav[3].data
-    lh_rot = left_handed.inav[3].data
-    mv_rot = multi_vector.inav[3].data
-
-    np.testing.assert_almost_equal(xy_rot, rh_rot, decimal=2)  # rotations
-    np.testing.assert_almost_equal(xy_rot, lh_rot, decimal=2)  # rotations
-    np.testing.assert_almost_equal(xy_rot, mv_rot, decimal=2)  # rotations
-
-
-@pytest.mark.skip(reason="basis change functionality not yet implemented")
-def test_trace(xy_vectors, right_handed, multi_vector):  # pragma: no cover
-    """
-    Basis does effect strain measurement, but we can simply calculate suitable invarients.
-    See https://en.wikipedia.org/wiki/Infinitesimal_strain_theory for details.
-    """
-    np.testing.assert_almost_equal(
-        np.add(
-            xy_vectors.inav[0].data, xy_vectors.inav[1].data), np.add(
-            right_handed.inav[0].data, right_handed.inav[1].data), decimal=2)
-    np.testing.assert_almost_equal(
-        np.add(
-            xy_vectors.inav[0].data, xy_vectors.inav[1].data), np.add(
-            multi_vector.inav[0].data, multi_vector.inav[1].data), decimal=2)

pyxem/tests/test_signals/test_strain_map.py

@@ -0,0 +1,119 @@
+# -*- coding: utf-8 -*-
+# Copyright 2017-2019 The pyXem developers
+#
+# This file is part of pyXem.
+#
+# pyXem is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# pyXem is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with pyXem.  If not, see <http://www.gnu.org/licenses/>.
+
+import hyperspy.api as hs
+import pytest
+import numpy as np
+from pyxem.tests.test_generators.test_displacement_gradient_tensor_generator import generate_test_vectors
+from pyxem.generators.displacement_gradient_tensor_generator import get_DisplacementGradientMap
+from pyxem.signals.strain_map import StrainMap, _get_rotation_matrix
+
+
+@pytest.fixture()
+def Displacement_Grad_Map():
+    xy = np.asarray([[1, 0], [0, 1]])
+    deformed = hs.signals.Signal2D(generate_test_vectors(xy))
+    D = get_DisplacementGradientMap(deformed, xy)
+    return D
+
+
+def test_rotation_matrix_formation():
+    x_new = [np.random.rand(), np.random.rand()]
+    R = _get_rotation_matrix(x_new)
+    ratio_array = np.divide(x_new, np.matmul(R, [1, 0]))
+    assert np.allclose(ratio_array[0], ratio_array[1])
+
+
+def test__init__(Displacement_Grad_Map):
+    strain_map = Displacement_Grad_Map.get_strain_maps()
+    assert strain_map.axes_manager.navigation_size == 4
+
+def test_signal_axes_carry_through(Displacement_Grad_Map):
+    """ A strain map that is calibrated, should stay calibrated when we change basis """
+    strain_map = Displacement_Grad_Map.get_strain_maps()
+    strain_map.axes_manager.signal_axes[1].units = 'nm'
+    strain_map.axes_manager.signal_axes[0].scale = 19
+    strain_alpha = strain_map.rotate_strain_basis([np.random.rand(), np.random.rand()])
+    assert strain_alpha.axes_manager.signal_axes[1].units == 'nm'
+    assert strain_alpha.axes_manager.signal_axes[0].scale == 19
+
+
+""" These are change of basis tests """
+
+
+def test_something_changes(Displacement_Grad_Map):
+    oneone_strain_original = Displacement_Grad_Map.get_strain_maps()
+    local_D = Displacement_Grad_Map
+    strain_alpha = local_D.get_strain_maps()
+    oneone_strain_alpha = strain_alpha.rotate_strain_basis([np.random.rand(), np.random.rand()])
+    assert not np.allclose(oneone_strain_original.data, oneone_strain_alpha.data, atol=0.01)
+
+
+def test_90_degree_rotation(Displacement_Grad_Map):
+    oneone_strain_original = Displacement_Grad_Map.get_strain_maps()
+    local_D = Displacement_Grad_Map
+    strain_alpha = local_D.get_strain_maps()
+    oneone_strain_alpha = strain_alpha.rotate_strain_basis([0, 1])
+    assert np.allclose(oneone_strain_original.inav[2:].data, oneone_strain_alpha.inav[2:].data, atol=0.01)
+    assert np.allclose(oneone_strain_original.inav[0].data, oneone_strain_alpha.inav[1].data, atol=0.01)
+    assert np.allclose(oneone_strain_original.inav[1].data, oneone_strain_alpha.inav[0].data, atol=0.01)
+
+
+def test_going_back_and_forward_between_bases(Displacement_Grad_Map):
+    """ Checks that going via an intermediate strain map doesn't give incorrect answers"""
+    strain_original = Displacement_Grad_Map.get_strain_maps()
+    local_D = Displacement_Grad_Map
+    temp_strain = local_D.get_strain_maps()
+    temp_strain = temp_strain.rotate_strain_basis([np.random.rand(), np.random.rand()])
+    fixed_xnew = [3.1, 4.1]
+    alpha = strain_original.rotate_strain_basis(fixed_xnew)
+    beta = temp_strain.rotate_strain_basis(fixed_xnew)
+    assert np.allclose(alpha, beta, atol=0.01)
+
+
+def test_rotation(Displacement_Grad_Map):
+    """
+    We should always measure the same rotations, regardless of basis
+    """
+    local_D = Displacement_Grad_Map
+    original = local_D.get_strain_maps()
+    rotation_alpha = original.rotate_strain_basis([np.random.rand(), np.random.rand()])
+    rotation_beta = original.rotate_strain_basis([np.random.rand(), -np.random.rand()])
+
+    # check the functionality has left invarient quantities invarient
+    np.testing.assert_almost_equal(original.inav[3].data, rotation_alpha.inav[3].data, decimal=2)  # rotations
+    np.testing.assert_almost_equal(original.inav[3].data, rotation_beta.inav[3].data, decimal=2)  # rotations
+
+
+def test_trace(Displacement_Grad_Map):
+    """
+    Basis does effect strain measurement, but we can simply calculate suitable invariants.
+    See https://en.wikipedia.org/wiki/Infinitesimal_strain_theory for details.
+    """
+
+    local_D = Displacement_Grad_Map
+    original = local_D.get_strain_maps()
+    rotation_alpha = original.rotate_strain_basis([1.3, +1.9])
+    rotation_beta = original.rotate_strain_basis([1.7, -0.3])
+
+    np.testing.assert_almost_equal(np.add(original.inav[0].data, original.inav[1].data),
+                                   np.add(rotation_alpha.inav[0].data, rotation_alpha.inav[1].data),
+                                   decimal=2)
+    np.testing.assert_almost_equal(np.add(original.inav[0].data, original.inav[1].data),
+                                   np.add(rotation_beta.inav[0].data, rotation_beta.inav[1].data),
+                                   decimal=2)