Merge 93f5d56 into accb829

diffsims/generators/rotation_list_generators.py

@@ -25,19 +25,75 @@
 from itertools import product
 
 from transforms3d.euler import euler2axangle, axangle2euler
+from transforms3d.euler import axangle2euler, euler2axangle, euler2mat
+from transforms3d.quaternions import quat2axangle, axangle2quat, mat2quat, qmult
+
+from diffsims.utils.rotation_conversion_utils import *
+from diffsims.utils.vector_utils import vectorised_spherical_polars_to_cartesians
+
+# Defines the maximum rotation angles [theta_max,psi_max,psi_min] associated with the
+# corners of the symmetry reduced region of the inverse pole figure for each crystal system.
+crystal_system_dictionary = {
+    "cubic": [45, 54.7, 0],
+    "hexagonal": [45, 90, 26.565],
+    "trigonal": [45, 90, -116.5],
+    "tetragonal": [45, 90, 0],
+    "orthorhombic": [90, 90, 0],
+    "monoclinic": [90, 0, -90],
+    "triclinic": [180, 360, 0],
+}
+
+def _rotate_axangle(Axangles, new_center):
+    """
+    Rotates a series of orientation described by axangle to a new center
+
+    Parameters
+    ----------
+    Axangles : diffsims.Axangles
+        Pre-rotation
+    new_center : (alpha,beta,gamma)
+        The location of the (0,0,0) rotation as an rzxz euler angle
+
+    Returns
+    -------
+    AxAngles : diffsims.Axangles
+        Rotated
+    See Also
+    --------
+    generators.get_local_grid
+    """
 
-from diffsims.utils.rotation_conversion_utils import Euler
-from diffsims.utils.fundamental_zone_utils import (
-    get_proper_point_group_string,
-    reduce_to_fundamental_zone,
-)
-from diffsims.utils.gridding_utils import (
-    create_linearly_spaced_array_in_rzxz,
-    rotate_axangle,
-    _create_advanced_linearly_spaced_array_in_rzxz,
-    get_beam_directions,
-    beam_directions_to_euler_angles,
-)
+    quats = Axangles.to_Quat()
+    q = mat2quat(rotation_matrix_from_euler_angles((new_center)))
+    stored_quats = vectorised_qmult(q, quats)
+
+    return AxAngle.from_Quat(stored_quats)
+
+
+def _beam_directions_to_euler_angles(points_in_cartesians):
+    """
+    Converts an array of cartesians (x,y,z unit basis vectors) to the euler angles that would take [0,0,1] to [x,y,z]
+    Parameters
+    ----------
+    points_in_cartesians :
+         Generally output from get_beam_directions()
+    Returns
+    -------
+    diffsims.Euler :
+         The appropriate euler angles
+    """
+    axes = np.cross(
+        [0, 0, 1], points_in_cartesians
+    )  # in unit cartesians so this is fine, [0,0,1] returns [0,0,0]
+    norms = np.linalg.norm(axes, axis=1).reshape(-1, 1)
+    angle = np.arcsin(norms)
+
+    normalised_axes = np.ones_like(axes)
+    np.divide(axes, norms, out=normalised_axes, where=norms != 0)
+
+    np_axangles = np.hstack((normalised_axes, angle.reshape((-1, 1))))
+    eulers = AxAngle(np_axangles).to_Euler(axis_convention="rzxz")
+    return eulers
 
 
 def _returnable_eulers_from_axangle(grid, axis_convention, round_to):
@@ -63,63 +119,14 @@ def get_fundamental_zone_grid(space_group_number, resolution):
     Returns
     -------
     rotation_list : list of tuples
-    """
-    zone_string = get_proper_point_group_string(space_group_number)
-    raw_grid = create_linearly_spaced_array_in_rzxz(
-        resolution
-    )  # see discussion in diffsims/#50
-    raw_grid_axangle = raw_grid.to_AxAngle()
-    fz_grid_axangle = reduce_to_fundamental_zone(raw_grid_axangle, zone_string)
-    return _returnable_eulers_from_axangle(fz_grid_axangle, "rzxz", round_to=2)
-
 
-def get_grid_stereographic(crystal_system, resolution, equal="angle"):
+    Notes
+    -----
     """
-    Creates a rotation list by determining the beam directions within the symmetry reduced
-    region of the inverse pole figure, corresponding to the specified crystal system, and
-    combining this with rotations about the beam direction at a given resolution.
-
-    Parameters
-    ----------
-    crytal_system : str
-        'cubic', 'hexagonal', 'trigonal', 'tetragonal', 'orthorhombic', 'monoclinic' and 'triclinic'
-
-    resolution : float
-        The maximum misorientation between rotations in the list, as defined according to
-        the parameter 'equal'. Specified as an angle in degrees.
-    equal : str
-        'angle' or 'area'. If 'angle', the misorientation is calculated between each beam direction
-        and its nearest neighbour(s). If 'area', the density of points is as in the equal angle case
-        but each point covers an equal area.
-
-    Returns
-    -------
-    rotation_list : list of tuples
-        List of rotations
-    """
-    beam_directions_rzxz = beam_directions_to_euler_angles(
-        get_beam_directions(crystal_system, resolution, equal=equal)
-    )
-    beam_directions_szxz = beam_directions_rzxz.to_AxAngle().to_Euler(
-        axis_convention="szxz"
-    )  # convert to high speed convention
-
-    # drop in all the inplane rotations to form z
-    alpha = beam_directions_szxz.data[:, 0]
-    beta = beam_directions_szxz.data[:, 1]
-    in_plane = np.arange(0, 360, resolution)
-
-    ipalpha = np.asarray(list(product(alpha, np.asarray(in_plane))))
-    ipbeta = np.asarray(list(product(beta, np.asarray(in_plane))))
-    z = np.hstack((ipalpha[:, 0].reshape((-1, 1)), ipbeta))
-
-    raw_grid = Euler(z, axis_convention="szxz")
-    grid_rzxz = raw_grid.to_AxAngle().to_Euler(
-        axis_convention="rzxz"
-    )  # convert back Bunge convention to return
-    rotation_list = grid_rzxz.to_rotation_list(round_to=2)
-    return rotation_list
-
+    #raise Deprecation warning
+    #g get grid from orix
+    #convert_to_rotation_list()
+    return None
 
 def get_local_grid(center, max_rotation, resolution):
     """
@@ -141,16 +148,18 @@ def get_local_grid(center, max_rotation, resolution):
     -------
     rotation_list : list of tuples
     """
-    raw_grid = _create_advanced_linearly_spaced_array_in_rzxz(
-        resolution, 360, max_rotation + 10, 360
-    )
-    raw_grid_axangle = raw_grid.to_AxAngle()
-    raw_grid_axangle.remove_large_rotations(np.deg2rad(max_rotation))
-    if not np.all(np.asarray(center) == 0):
-        raw_grid_axangle = rotate_axangle(raw_grid_axangle, center)
+    #raise Deprecation warning
+    #g get grid from orix
+    #convert_to_rotation_list()
+    return None
 
-    return _returnable_eulers_from_axangle(raw_grid_axangle, "rzxz", round_to=2)
 
+def get_grid_stereographic(crystal_system, resolution, equal="angle"):
+    """
+    This functionality is deprecated. The following outline is only given to
+    aid dev work
+    """
+    return get_fundamental_zone_grid(1,resolution)
 
 def get_grid_around_beam_direction(beam_rotation, resolution, angular_range=(0, 360)):
     """
@@ -199,3 +208,98 @@ def get_grid_around_beam_direction(beam_rotation, resolution, angular_range=(0,
     grid_rzxz = raw_grid.to_AxAngle().to_Euler(axis_convention="rzxz")
     rotation_list = grid_rzxz.to_rotation_list(round_to=2)
     return rotation_list
+
+
+#rename this to, get_beam_direction_grid()
+def get_beam_directions(crystal_system, resolution, equal="angle"):
+    """
+    Produces an array of beam directions, evenly (see equal argument) spaced that lie within the streographic
+    triangle of the relevant crystal system.
+
+    Parameters
+    ----------
+    crystal_system : str
+        Allowed are: 'cubic','hexagonal','trigonal','tetragonal','orthorhombic','monoclinic','triclinic'
+
+    resolution : float
+        An angle in degrees. If the 'equal' option is set to 'angle' this is the misorientation between a
+        beam direction and its nearest neighbour(s). For 'equal'=='area' the density of points is as in
+        the equal angle case but each point covers an equal area
+
+    equal : str
+        'angle' (default) or 'area'
+
+    Returns
+    -------
+    points_in_cartesians : np.array (N,3)
+        Rows are x,y,z where z is the 001 pole direction.
+    Notes
+    -----
+    For all cases: The input 'resolution' may differ slightly from the expected value. This is so that each of the corners
+    of the streographic triangle are included. Actual 'resolution' will always be equal to or higher than the input resolution. As
+    an example, if resolution is set to 4 to cover a range [0,90] we can't include both endpoints. The code makes 23 steps
+    of 3.91 degrees instead.
+
+    For the cubic case: Each edge of the streographic triangle will behave as expected. The region above the (1,0,1), (1,1,1) edge
+    will (for implementation reasons) be slightly more densly packed than the wider region.
+    """
+    theta_max, psi_max, psi_min = crystal_system_dictionary[crystal_system]
+
+    # see docstrings for np.arange, np.linspace has better endpoint handling
+    steps_theta = int(np.ceil((theta_max - 0) / resolution))
+    steps_psi = int(np.ceil((psi_max - psi_min) / resolution))
+    theta = np.linspace(
+        0, np.deg2rad(theta_max), num=steps_theta
+    )  # radians as we're about to make spherical polar cordinates
+    if equal == "area":
+        # http://mathworld.wolfram.com/SpherePointPicking.html
+        v_1 = (1 + np.cos(np.deg2rad(psi_max))) / 2
+        v_2 = (1 + np.cos(np.deg2rad(psi_min))) / 2
+        v_array = np.linspace(min(v_1, v_2), max(v_1, v_2), num=steps_psi)
+        psi = np.arccos(2 * v_array - 1)  # in radians
+    elif equal == "angle":
+        # now in radians as we're about to make spherical polar cordinates
+        psi = np.linspace(np.deg2rad(psi_min), np.deg2rad(psi_max), num=steps_psi)
+
+    psi_theta = np.asarray(list(product(psi, theta)))
+    r = np.ones((psi_theta.shape[0], 1))
+    points_in_spherical_polars = np.hstack((r, psi_theta))
+
+    # keep only theta ==0 psi ==0, do this with np.abs(theta) > 0 or psi == 0 - more generally use the smallest psi value
+    points_in_spherical_polars = points_in_spherical_polars[
+        np.logical_or(
+            np.abs(psi_theta[:, 1]) > 0, psi_theta[:, 0] == np.min(psi_theta[:, 0])
+        )
+    ]
+    points_in_cartesians = vectorised_spherical_polars_to_cartesians(
+        points_in_spherical_polars
+    )
+
+    if crystal_system == "cubic":
+        # add in the geodesic that runs [1,1,1] to [1,0,1]
+        v1 = np.divide([1, 1, 1], np.sqrt(3))
+        v2 = np.divide([1, 0, 1], np.sqrt(2))
+
+        def cubic_corner_geodesic(t):
+            # https://math.stackexchange.com/questions/1883904/a-time-parameterization-of-geodesics-on-the-sphere
+            w = v2 - np.multiply(np.dot(v1, v2), v1)
+            w = np.divide(w, np.linalg.norm(w))
+            # return in cartesians with t_end = np.arccos(np.dot(v1,v2))
+            return np.add(
+                np.multiply(np.cos(t.reshape(-1, 1)), v1),
+                np.multiply(np.sin(t.reshape(-1, 1)), w),
+            )
+
+        t_list = np.linspace(0, np.arccos(np.dot(v1, v2)), num=steps_theta)
+        geodesic = cubic_corner_geodesic(t_list)
+        points_in_cartesians = np.vstack((points_in_cartesians, geodesic))
+        # the great circle (from [1,1,1] to [1,0,1]) forms a plane (with the
+        # origin), points on the same side as (0,0,1) are safe, the others are not
+        plane_normal = np.cross(
+            v2, v1
+        )  # dotting this with (0,0,1) gives a positive number
+        points_in_cartesians = points_in_cartesians[
+            np.dot(points_in_cartesians, plane_normal) >= 0
+        ]  # 0 is the points on the geodesic
+
+    return points_in_cartesians

diffsims/libraries/structure_library.py

@@ -80,6 +80,7 @@ def from_orientation_lists(cls, identifiers, structures, orientations):
         """
         return cls(identifiers, structures, orientations)
 
+    # TODO: this should be from space groups
     @classmethod
     def from_crystal_systems(
         cls, identifiers, structures, systems, resolution, equal="angle"
@@ -114,7 +115,7 @@ def get_library_size(self, to_print=False):
         Returns the the total number of orientations in the
         current StructureLibrary object. Will also print the number of orientations
         for each identifier in the library if the to_print==True
-        
+
         Parameters
         ----------
         to_print : bool

diffsims/tests/test_generators/test_rotation_list_generator.py

@@ -21,64 +21,16 @@
 from diffsims.generators.rotation_list_generators import (
     get_local_grid,
     get_grid_around_beam_direction,
-    get_fundamental_zone_grid,
-    get_grid_stereographic,
-)
-from diffsims.utils.rotation_conversion_utils import Euler
-
-
-@pytest.mark.parametrize("center", [(0.0, 0.0, 0.0), (0.0, 10.0, 0.0)])
-def test_get_local_grid(center):
-    grid = get_local_grid(center, 10, 2)
-    assert isinstance(grid, list)
-    assert isinstance(grid[0], tuple)
-    assert center in grid
-    center_plus_2 = (center[0], center[1], center[2] + 2)
-    assert center_plus_2 in grid
+    get_fundamental_zone_grid
+    )
 
+@pytest.mark.skip(reason="Not implemented")
+def test_calls_to_orix():
+    _ = get_local_grid((0,0,0),31,15)
+    _ = get_fundamental_zone_grid(space_group_number, resolution=20)
 
 def test_get_grid_around_beam_direction():
     grid_simple = get_grid_around_beam_direction([1, 1, 1], 1, (0, 360))
     assert isinstance(grid_simple, list)
     assert isinstance(grid_simple[0], tuple)
     assert len(grid_simple) == 360
-
-
-@pytest.mark.parametrize("space_group_number", [1, 3, 30, 190, 215, 229])
-def test_get_fundamental_zone_grid(space_group_number):
-    grid_narrow = get_fundamental_zone_grid(space_group_number, resolution=4)
-    grid_wide = get_fundamental_zone_grid(space_group_number, resolution=8)
-    assert (len(grid_narrow) / len(grid_wide)) > (2 ** 3) - 1  # lower bounds
-    assert (len(grid_narrow) / len(grid_wide)) < (2 ** 3) + 1  # upper bounds
-
-
-def test_get_grid_stereographic():
-    grid = get_grid_stereographic("tetragonal", 1, equal="angle")
-    assert (0, 0, 0) in grid
-    grid_four_times_as_many = get_grid_stereographic("orthorhombic", 1, equal="angle")
-
-    # for equal angle you wouldn't expect perfect ratios
-    assert len(grid_four_times_as_many) / len(grid) > 1.9
-    assert len(grid_four_times_as_many) / len(grid) < 2.1
-
-
-@pytest.mark.skip(reason="This tests a theoretical underpinning of the code")
-def test_small_angle_shortcut():  # pragma: no cover
-    """ Demonstrates that cutting larger 'out of plane' in euler space doesn't
-    effect the result """
-
-    def process_angles(raw_angles, max_rotation):
-        raw_angles = raw_angles.to_AxAngle()
-        raw_angles.remove_large_rotations(np.deg2rad(max_rotation))
-        return raw_angles
-
-    max_rotation = 20
-    lsa = create_linearly_spaced_array_in_rzxz(2)
-    alsa = _create_advanced_linearly_spaced_array_in_rzxz(
-        2, 360, max_rotation + 10, 360
-    )
-
-    long_true_way = process_angles(lsa, max_rotation)
-    quick_way = process_angles(alsa, max_rotation)
-
-    assert long_true_way.data.shape == quick_way.data.shape