Lambert projection

doc/changelog.rst

@@ -17,11 +17,15 @@ Unreleased
 Contributors
 ------------
 - Håkon Wiik Ånes
+- Lars Andreas Hastad Lervik
 
 Added
 -----
 - Reader for EMsoft's simulated EBSD patterns returned by their ``EMEBSD.f90``
   program. (`#202 <https://github.com/pyxem/kikuchipy/pull/202>`_)
+- Modified Lambert mapping, and its inverse, from points on the unit sphere to a
+  2D square grid, as implemented in Callahan and De Graef (2013).
+  (`#214 <https://github.com/pyxem/kikuchipy/pull/214>`_)
 
 Changed
 -------

kikuchipy/projections/__init__.py

@@ -20,12 +20,14 @@
     ebsd_projections,
     gnomonic_projection,
     hesse_normal_form,
+    lambert_projection,
     spherical_projection,
 )
 
 __all__ = [
     "ebsd_projections",
-    "gnomonic_projection.py",
+    "gnomonic_projection",
     "hesse_normal_form",
-    "spherical_projection.py",
+    "lambert_projection",
+    "spherical_projection",
 ]

kikuchipy/projections/gnomonic_projection.py

@@ -27,7 +27,9 @@
 class GnomonicProjection(SphericalProjection):
     """Gnomonic projection of a vector as implemented in MTEX."""
 
+    @classmethod
     def project(self, v: Union[Vector3d, np.ndarray]) -> np.ndarray:
+        """Convert from Cartesian to the Gnomonic projection."""
         polar = super().project(v)
         theta, phi = polar[..., 0], polar[..., 1]
 
@@ -53,6 +55,7 @@ def project(self, v: Union[Vector3d, np.ndarray]) -> np.ndarray:
 
     @staticmethod
     def iproject(x: np.ndarray, y: np.ndarray) -> Vector3d:
+        """Convert from the Gnomonic projection to Cartesian coordinates."""
         theta = np.arctan(np.sqrt(x ** 2 + y ** 2))
-        phi = np.atan2(y, x)
+        phi = np.arctan2(y, x)
         return Vector3d.from_polar(theta=theta, phi=phi)

kikuchipy/projections/lambert_projection.py

@@ -0,0 +1,127 @@
+# -*- coding: utf-8 -*-
+# Copyright 2019-2020 The kikuchipy developers
+#
+# This file is part of kikuchipy.
+#
+# kikuchipy 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.
+#
+# kikuchipy 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 kikuchipy. If not, see <http://www.gnu.org/licenses/>.
+
+from typing import Union
+
+import numpy as np
+from orix.vector import Vector3d
+
+from kikuchipy.projections.gnomonic_projection import GnomonicProjection
+
+
+class LambertProjection:
+    """Lambert projection of a vector [Callahan2013]_."""
+
+    @classmethod
+    def project(cls, v: Union[Vector3d, np.ndarray]) -> np.ndarray:
+        """Convert (n, 3) vector from Cartesian to the Lambert projection."""
+        if isinstance(v, Vector3d):
+            w = v.unit.data
+            x = w[..., 0]
+            y = w[..., 1]
+            z = w[..., 2]
+        else:
+            norm = np.sqrt(
+                np.sum(np.square([v[..., 0], v[..., 1], v[..., 2]]), axis=0)
+            )
+            x = v[..., 0] / norm
+            y = v[..., 1] / norm
+            z = v[..., 2] / norm
+
+        # Arrays used in both setting X and Y
+        sqrt_z = np.sqrt(2 * (1 - z))
+        sign_x = np.sign(x)
+        sign_y = np.sign(y)
+        abs_yx = abs(y) <= abs(x)
+
+        # Reusable constants
+        sqrt_pi = np.sqrt(np.pi)
+        sqrt_pi_half = sqrt_pi / 2
+        two_over_sqrt_pi = 2 / sqrt_pi
+
+        # Equations 10a and 10b from Callahan and De Graef (2013)
+        X = np.where(
+            abs_yx,
+            sign_x * sqrt_z * sqrt_pi_half,
+            sign_y * sqrt_z * (two_over_sqrt_pi * np.arctan2(x, y)),
+        )
+        Y = np.where(
+            abs_yx,
+            sign_x * sqrt_z * (two_over_sqrt_pi * np.arctan2(y, x)),
+            sign_y * sqrt_z * sqrt_pi_half,
+        )
+
+        return np.column_stack((X, Y))
+
+    @staticmethod
+    def iproject(xy: np.ndarray) -> Vector3d:
+        """Convert (n, 2) array from Lambert to Cartesian coordinates."""
+        X = xy[..., 0]
+        Y = xy[..., 1]
+
+        # Arrays used in setting x and y
+        true_term = Y * np.pi / (4 * X)
+        false_term = X * np.pi / (4 * Y)
+        abs_yx = abs(Y) <= abs(X)
+        c_x = _eq_c(X)
+        c_y = _eq_c(Y)
+
+        # Equations 8a and 8b from Callahan and De Graef (2013)
+        x = np.where(
+            abs_yx,
+            c_x * np.cos(true_term),
+            c_y * np.sin(false_term),
+        )
+        y = np.where(
+            abs_yx,
+            c_x * np.sin(true_term),
+            c_y * np.cos(false_term),
+        )
+        z = np.where(
+            abs_yx,
+            1 - (2 * (X ** 2)) / np.pi,
+            1 - (2 * (Y ** 2)) / np.pi,
+        )
+
+        return Vector3d(np.column_stack((x, y, z)))
+
+    @staticmethod
+    def lambert_to_gnomonic(xy: np.ndarray) -> np.ndarray:
+        """Convert (n,2) array from Lambert via Cartesian coordinates to
+        Gnomonic."""
+        # These two functions could probably be combined into 1 to decrease
+        # runtime
+        vec = LambertProjection.iproject(xy)
+        xy = GnomonicProjection.project(vec)
+        return xy
+
+    @staticmethod
+    def gnomonic_to_lambert(xy: np.ndarray) -> np.ndarray:
+        """Convert (n,2) array from Gnomonic via Cartesian coordinates to
+        Lambert."""
+        # These two functions could probably be combined into 1 to decrease
+        # runtime
+        vec = GnomonicProjection.iproject(xy[..., 0], xy[..., 1])
+        xy = LambertProjection.project(vec)
+        return xy
+
+
+def _eq_c(p: Union[np.ndarray, float, int]) -> Union[np.ndarray, float, int]:
+    """Private function used inside LambertProjection.iproject to increase
+    readability."""
+    return (2 * p / np.pi) * np.sqrt(np.pi - p ** 2)

kikuchipy/projections/spherical_projection.py

@@ -28,6 +28,7 @@ class SphericalProjection:
 
     spherical_region = SphericalRegion([0, 0, 1])
 
+    @classmethod
     def project(self, v: Union[Vector3d, np.ndarray]) -> np.ndarray:
         """Convert from cartesian to spherical coordinates."""
         polar = get_polar(v)

kikuchipy/projections/tests/__init__.py

@@ -0,0 +1,17 @@
+# -*- coding: utf-8 -*-
+# Copyright 2019-2020 The kikuchipy developers
+#
+# This file is part of kikuchipy.
+#
+# kikuchipy 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.
+#
+# kikuchipy 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 kikuchipy. If not, see <http://www.gnu.org/licenses/>.

kikuchipy/projections/tests/test_lambert_projection.py

@@ -0,0 +1,92 @@
+# -*- coding: utf-8 -*-
+# Copyright 2019-2020 The kikuchipy developers
+#
+# This file is part of kikuchipy.
+#
+# kikuchipy 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.
+#
+# kikuchipy 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 kikuchipy. If not, see <http://www.gnu.org/licenses/>.
+
+import numpy as np
+from orix.vector import Vector3d
+import pytest
+
+from kikuchipy.projections.lambert_projection import LambertProjection, _eq_c
+
+
+class TestLambertProjection:
+    def test_project_vector3d(self):
+        """Works for Vector3d objects with single and multiple vectors"""
+
+        vector_one = Vector3d((0.578, 0.578, 0.578))
+        output_a = LambertProjection.project(vector_one)
+        expected_a = np.array((0.81417, 0.81417))
+
+        assert output_a[..., 0][0] == pytest.approx(expected_a[..., 0], rel=1e4)
+        assert output_a[..., 1][0] == pytest.approx(expected_a[..., 1], rel=1e4)
+
+        vector_two = Vector3d(
+            np.array(
+                [[0.578, 0.578, 0.578], [0, 0.707, 0.707], [0.707, 0, 0.707]]
+            )
+        )
+        output_b = LambertProjection.project(vector_two)
+
+        expected_b = np.array(
+            (
+                (0.814172, 0.814172, 0.814172),
+                (0, 0, 0),
+                (0.6784123, 0, 0.6784123),
+            )
+        )
+
+        assert output_a[..., 0][0] == pytest.approx(expected_a[..., 0], rel=1e4)
+        assert output_a[..., 1][0] == pytest.approx(expected_a[..., 1], rel=1e4)
+
+    def test_project_ndarray(self):
+        "Works for numpy ndarrays"
+        ipt = np.array((0.578, 0.578, 0.578))
+        output = LambertProjection.project(ipt)
+        expected = np.array((0.81417, 0.81417))
+        assert output[..., 0][0] == pytest.approx(expected[..., 0], rel=1e4)
+
+    def test_iproject(self):
+        """Conversion from Lambert to Cartesian coordinates works"""
+        vec = np.array((0.81417, 0.81417))
+        expected = Vector3d((0.57702409, 0.577, 0.578))
+        output = LambertProjection.iproject(vec)
+        assert np.allclose(output.x.data[0], expected.x.data[0], rtol=1e-05)
+
+    def test_eq_c(self):
+        """Helper function works"""
+        arr_1 = np.array((0.578, 0.578, 0.578))
+        value = arr_1[..., 0]
+        output_arr = np.array((0.61655, 0.61655, 0.61655))
+        expected = output_arr[..., 0]
+        output = _eq_c(value)
+        assert output == pytest.approx(expected, rel=1e-5)
+
+    def test_lambert_to_gnomonic(self):
+        """Conversion from Lambert to Gnomonic works"""
+        vec = np.array(
+            (0.81417, 0.81417)
+        )  # Should give x,y,z = 1/sqrt(3) (1, 1, 1)
+        output = LambertProjection.lambert_to_gnomonic(vec)
+        expected = np.array((1, 1))
+        assert output[..., 0][0] == pytest.approx(expected[..., 0], abs=1e-2)
+
+    def test_gnomonic_to_lambert(self):
+        """Conversion from Gnomonic to Lambert works"""
+        vec = np.array((1, 1))  # Similar to the case above
+        output = LambertProjection.gnomonic_to_lambert(vec)
+        expected = np.array((0.81417, 0.81417))
+        assert output[..., 0][0] == pytest.approx(expected[..., 0], rel=1e-3)

kikuchipy/release.py

@@ -18,7 +18,11 @@
 
 author = "kikuchipy developers"
 copyright = "Copyright 2019-2020, kikuchipy"
-credits = ["Håkon Wiik Ånes", "Tina Bergh"]
+credits = [
+    "Håkon Wiik Ånes",
+    "Tina Bergh",
+    "Lars Andreas Hastad Lervik",
+]
 license = "GPLv3+"
 maintainer = "Håkon Wiik Ånes"
 maintainer_email = "hakon.w.anes@ntnu.no"

kikuchipy/signals/ebsd.py

@@ -1531,3 +1531,4 @@ def get_decomposition_model_write(
         # Delete temporary files
         os.remove(file_learn)
         gc.collect()  # Don't sink
+