Rotation lists from streographic gridding

diffsims/generators/rotation_list_generators.py

@@ -22,16 +22,18 @@
 
 import numpy as np
 import warnings
+from itertools import product
 
 from diffsims.utils.rotation_conversion_utils import Euler
 from diffsims.utils.fundemental_zone_utils import get_proper_point_group_string, reduce_to_fundemental_zone
 from diffsims.utils.gridding_utils import create_linearly_spaced_array_in_rzxz,rotate_axangle, \
                                           _create_advanced_linearly_spaced_array_in_rzxz, \
-                                          _get_rotation_to_beam_direction
+                                          _get_rotation_to_beam_direction, get_beam_directions, beam_directions_to_euler_angles
 
 
 def _returnable_eulers_from_axangle(grid,axis_convention,round_to):
-    """ Converts a grid of orientations in axis-angle space to Euler angles following a user specified convention and rounding."""
+    """ Converts a grid of orientations in axis-angle space to Euler
+    angles following a user specified convention and rounding."""
     eulers = grid.to_Euler(axis_convention=axis_convention)
     rotation_list = eulers.to_rotation_list(round_to=round_to)
     return rotation_list
@@ -58,14 +60,55 @@ def get_fundemental_zone_grid(space_group_number, resolution):
     fz_grid_axangle = reduce_to_fundemental_zone(raw_grid_axangle, zone_string)
     return _returnable_eulers_from_axangle(fz_grid_axangle,'rzxz',round_to=2)
 
+def get_grid_streographic(crystal_system,resolution,equal='angle'):
+    """
+    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
+
 
 def get_local_grid(center, max_rotation, resolution):
     """
     Creates a rotation list for the rotations within max_rotation of center at a given rotation.
 
     Parameters
     ----------
-    center : 3 angle tuple
+    center : tuple
         The orientation that acts as the center of the grid, as euler angles specified in the
         'rzxz' convention (degrees)
 
@@ -103,15 +146,14 @@ def get_grid_around_beam_direction(beam_direction,resolution, angular_range=(0,
     angular_range : tuple
         The minimum (included) and maximum (excluded) rotation around the beam direction to be included
 
-    cubic : bool (Default=False)
+    cubic : bool
         This only works for cubic systems at the present, when False this raises a warning, set to
-        True to supress said warning
+        True to supress said warning. The default is False
 
     Returns
     -------
     rotation_list : list of tuples
     """
-    from itertools import product
 
     if not cubic:
         warnings.warn("This code only works for cubic systems at present")

diffsims/tests/test_generators/test_rotation_list_generator.py

@@ -18,7 +18,7 @@
 
 import pytest
 import numpy as np
-from diffsims.generators.rotation_list_generators import get_local_grid, get_grid_around_beam_direction,get_fundemental_zone_grid
+from diffsims.generators.rotation_list_generators import get_local_grid, get_grid_around_beam_direction, get_fundemental_zone_grid, get_grid_streographic
 from diffsims.utils.rotation_conversion_utils import Euler
 
 
@@ -41,6 +41,14 @@ def test_get_grid_around_beam_direction():
 def test_get_fundemental_zone_grid(space_group_number):
     grid = get_fundemental_zone_grid(space_group_number,resolution=3)
 
+def test_get_grid_streographic():
+    grid = get_grid_streographic('tetragonal',1,equal='angle')
+    assert (0,0,0) in grid
+    grid_four_times_as_many = get_grid_streographic('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

diffsims/tests/test_utils/test_gridding_utils.py

@@ -22,7 +22,7 @@
 from diffsims.utils.rotation_conversion_utils import AxAngle, Euler
 from diffsims.utils.fundemental_zone_utils import get_proper_point_group_string, reduce_to_fundemental_zone, numpy_bounding_plane
 from diffsims.utils.gridding_utils import create_linearly_spaced_array_in_rzxz, vectorised_qmult, \
-                                          _create_advanced_linearly_spaced_array_in_rzxz
+                                          _create_advanced_linearly_spaced_array_in_rzxz, get_beam_directions
 
 from transforms3d.quaternions import qmult
 
@@ -64,3 +64,48 @@ def test_select_fundemental_zone():
 def test_edge_case_numpy_bounding_plane():
     z = np.asarray([1,1,1,np.inf])
     numpy_bounding_plane(data=z,vector=[1,1,1],distance=1)
+
+""" This tests get_beam_directions """
+@pytest.mark.parametrize("crystal_system,expected_corners",
+                         [
+                         ['monoclinic',[(0,0,1),(0,1,0),(0,-1,0)]],
+                         ['orthorhombic',[(0,0,1),(1,0,0),(0,1,0)]],
+                         ['tetragonal',[(0,0,1),(1,0,0),(1,1,0)]],
+                         ['cubic',[(0,0,1),(1,0,1),(1,1,1)]],
+                         ['hexagonal',[(0,0,1),(2,1,0),(1,1,0)]],
+                         ['trigonal',[(0,0,1),(-2,-1,0),(1,1,0)]]
+                         ])
+def test_get_beam_directions_equal_angle(crystal_system,expected_corners):
+    z = get_beam_directions(crystal_system,1,'angle')
+    assert np.allclose(np.linalg.norm(z,axis=1),1)
+    for corner in expected_corners:
+        norm_corner = np.divide(corner,np.linalg.norm(corner))
+        assert np.any(np.isin(z,norm_corner))
+
+@pytest.mark.parametrize("crystal_system,expected_corners",
+                         [
+                         ['monoclinic',[(0,0,1),(0,1,0),(0,-1,0)]],
+                         ['orthorhombic',[(0,0,1),(1,0,0),(0,1,0)]],
+                         ['tetragonal',[(0,0,1),(1,0,0),(1,1,0)]],
+                         ['cubic',[(0,0,1),(1,0,1),(1,1,1)]],
+                         ['hexagonal',[(0,0,1),(2,1,0),(1,1,0)]],
+                         ['trigonal',[(0,0,1),(-2,-1,0),(1,1,0)]]
+                         ])
+def test_get_beam_directions_equal_area(crystal_system,expected_corners):
+    z = get_beam_directions(crystal_system,1,equal='area')
+    assert np.allclose(np.linalg.norm(z,axis=1),1)
+    for corner in expected_corners:
+        norm_corner = np.divide(corner,np.linalg.norm(corner))
+        assert np.any(np.isin(z,norm_corner))
+
+@pytest.mark.parametrize("crystal_system",['cubic','hexagonal'])
+def test_equal_area_same_as_equal_angle(crystal_system):
+    z_angle = get_beam_directions(crystal_system,1,equal='angle')
+    z_area = get_beam_directions(crystal_system,1,equal='area')
+    assert np.all(z_angle.shape==z_area.shape)
+
+def test_beam_directions_cubic():
+    # Following "Orientation precision of TEM-based orientation mapping techniques" - Morawiec et al, Ultramicroscopy 136,2014
+    z = get_beam_directions('cubic',1.6)
+    assert z.shape[0] > 950
+    assert z.shape[0] < 1050

diffsims/utils/gridding_utils.py

@@ -25,6 +25,17 @@
 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':[360,180,0]}
 
 
 def vectorised_qmult(q1, qdata):
@@ -94,15 +105,15 @@ def rotate_axangle(Axangles, new_center):
 
     Parameters
     ----------
-    Axangles :
-        Axangles in the correct class
+    Axangles : diffsims.Axangles
+        Pre-rotation
     new_center : (alpha,beta,gamma)
         The location of the (0,0,0) rotation as an rzxz euler angle
 
-    Returns 
+    Returns
     -------
-    AxAngles :
-
+    AxAngles : diffsims.Axangles
+        Rotated
     See Also
     --------
     generators.get_local_grid
@@ -160,7 +171,9 @@ def _create_advanced_linearly_spaced_array_in_rzxz(resolution, max_alpha, max_be
     diffsims.Euler
 
     """
-    steps_alpha = int(np.ceil((max_alpha - 0)/resolution)) #see docstrings for np.arange, np.linspace has better endpoint handling
+    # We use np.linspace rather than np.arange to get list of evenly spaced Euler
+    # angles due to better end point handling. Therefore convert "step_size" to a "num"
+    steps_alpha = int(np.ceil((max_alpha - 0)/resolution))
     steps_beta  = int(np.ceil((max_beta  - 0)/resolution))
     steps_gamma = int(np.ceil((max_gamma - 0)/resolution))
 
@@ -169,3 +182,98 @@ def _create_advanced_linearly_spaced_array_in_rzxz(resolution, max_alpha, max_be
     gamma = np.linspace(0, max_gamma, num=steps_gamma, endpoint=False)
     z = np.asarray(list(product(alpha, beta, gamma)))
     return Euler(z, axis_convention='rzxz')
+
+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]
+
+    steps_theta = int(np.ceil((theta_max - 0)/resolution)) #see docstrings for np.arange, np.linspace has better endpoint handling
+    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
+
+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]
+    angle = np.arcsin(np.linalg.norm(axes,axis=1))
+    normalised_axes = np.where(angle.reshape(-1,1) > 0, np.divide(axes,np.linalg.norm(axes,axis=1).reshape(-1,1)), axes)
+    np_axangles = np.hstack((normalised_axes,angle.reshape((-1,1))))
+    eulers = AxAngle(np_axangles).to_Euler(axis_convention='rzxz')
+    return eulers

diffsims/utils/vector_utils.py

@@ -65,3 +65,28 @@ def get_angle_cartesian(a, b):
     if denom == 0:
         return 0.0
     return math.acos(max(-1.0, min(1.0, np.dot(a, b) / denom)))
+
+def vectorised_spherical_polars_to_cartesians(z):
+    """
+    Converts an array of spherical polars into an array of
+    (x,y,z) = r(cos(psi)sin(theta),sin(psi)sin(theta),cos(theta))
+
+    Parameters
+    ----------
+    z : np.array
+        With rows of
+        r : the radius value, r = sqrt(x**2+y**2+z**2)
+        psi : The azimuthal angle generally (0,2pi])
+        theta : The elevation angle generally (0,pi)
+
+    Returns
+    -------
+    xyz : np.array
+        With rows of
+        x,y,z
+    """
+    r, psi, theta = z[:,0],z[:,1],z[:,2]
+    x = r * np.cos(psi) * np.sin(theta)
+    y = r * np.sin(psi) * np.sin(theta)
+    z = r * np.cos(theta)
+    return np.asarray([x,y,z]).T