Lambert projection (#214)

* Added dummy comment * Remove dummy comment * Added dummy comment * Remove dummy comment * Created lambert_projection.py for lambert projection * Updated __init__.py in projections * Created LambertProjection class in the same style as other projection styles with empty class methods. * Created TODOs on the empty methods * "to lambert" equations should now be in place. Brain needed to decide which path to take still TODO. * Fixed output vector for iproject method in Lambert projection * Created empty methods for the conversion between lambert and gnomonic projections * Implemented lambert_to_gnomonic method. Required changing GnomonicProjection.project to a class method. Might not be ok as is!! * Implemented gnomonic_to_lambert method. Required changing LambertProjection.project to a class method. Might not be ok as is!! * Fixed LambertProjection.project method to check the conditional and the method looks completed (not tested) * Fixed LambertProjection.iproject method to check the conditional and the method looks completed (not tested) * Fixed typos * Added explanation to eq_c method * Made eq_c a private method * Added Edge case fix for iproject method and optimize TOOD. {skip CI} * Moved/Deleted some inline comments. Added some TODOs {skip CI} * Added docstrings to every method in lambert_projection. {skip CI} * Added docstrings to every method in gnomonic_projection. {skip CI} * Started working briefly on some tests for the lambert projection module. {skip CI} * Added TODOs for critical methods missing (found while writing tests). {skip CI} * [skip ci] Fixed x, y, z assignments in LambertProjection.project as mentioned in #issuecomment-690713128 * [skip ci] LambertProjection.project now raises an error if it gets a vector that is not on the unit sphere. * [skip ci] LambertProjection.project now only uses numpy arrays. Pretty sure the for-loop is here to stay.. * [skip ci] LambertProjection.iproject now only uses numpy arrays. Pretty sure the for-loop is here to stay.. * [skip ci] LambertProjection.iproject now only uses numpy arrays. Pretty sure the for-loop is here to stay.. The method now also raises a ValueError if the X/Y input is greater than L. This is to keep it in sync with the Unit sphere of the project method. * [skip ci] Added TODO to assert that input vector in GnomicProjection.project is actually on the unit sphere. * [skip ci] Made the "on-sphere" checker less strict * [skip ci] Moved index variable to where it belongs :) * Made tests for every method. I expect them to fail initially, mostly due to floats * [skip ci] Fixed a typo in gnomonic_projection. np.atan2 -> np.arctan2 * [skip ci] Made SphericalProjection.project a classmethod * Made tests for every method, and fixed the code where it was needed. * [skip ci] Changed LambertProjection to only take one vector instead of the elements of the vector to make it consistent with the rest. Similar changes should probably be made in gnomonic_projection.py in the future. * [skip ci] Helper function skeletons made * [skip ci] Finished helper functions, and removed duplicated code * [skip ci] Finished writing tests * [skip ci] Updated type hints * Replaced for-loops with numpy.where() Updated methods to follow PEP8 Swapped some test methods with numpy.allclose(). Pytest.approx() still used mostly due to simplicity. [skip ci] * Updated docstring to fit [skip ci] * Deleted a comment, woke travis * Add @friedkitteh to credits, update changelog [skip ci] Signed-off-by: Håkon Wiik Ånes <hwaanes@gmail.com> * Updated docstrings [skip ci] * Fixed import statement in test_lambert_projection.py [skip ci] * Re-use arrays in LambertProjection.iproject [skip ci] * Updated docstring [skip ci] * Added comments [skip ci] * Updated docstrings [skip ci] * Updated docstrings [skip ci] * Optimized project and iproject methods [skip ci] * Optimized iproject method [skip ci] * Optimized iproject method [skip ci] * Ran "black ." to create a change to create a commit to ping travis Signed-off-by: Håkon Wiik Ånes <hwaanes@gmail.com> Co-authored-by: Håkon Wiik Ånes <hwaanes@gmail.com>
pyxem · Sep 15, 2020 · 80b98a9 · 80b98a9
1 parent b0846ef
commit 80b98a9
Show file tree

Hide file tree

Showing 8 changed files with 244 additions and 4 deletions.
diff --git a/doc/changelog.rst b/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
 -------

diff --git a/kikuchipy/projections/__init__.py b/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",
 ]
diff --git a/kikuchipy/projections/gnomonic_projection.py b/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)
diff --git a/kikuchipy/projections/lambert_projection.py b/kikuchipy/projections/lambert_projection.py
@@ -0,0 +1,117 @@
+# -*- 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)
diff --git a/kikuchipy/projections/spherical_projection.py b/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)

diff --git a/kikuchipy/projections/tests/__init__.py b/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/>.
diff --git a/kikuchipy/projections/tests/test_lambert_projection.py b/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)
diff --git a/kikuchipy/release.py b/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"