Merge pull request #104 from pc494/remove_gridding

Remove gridding

diffsims/generators/rotation_list_generators.py

@@ -24,133 +24,97 @@
 import warnings
 from itertools import product
 
-from transforms3d.euler import euler2axangle, axangle2euler
-
-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,
-)
-
-
-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."""
-    eulers = grid.to_Euler(axis_convention=axis_convention)
-    rotation_list = eulers.to_rotation_list(round_to=round_to)
-    return rotation_list
+from orix.sampling.sample_generators import get_sample_fundamental, get_sample_local
 
-
-def get_fundamental_zone_grid(space_group_number, resolution):
+from transforms3d.euler import euler2axangle, axangle2euler
+from transforms3d.euler import axangle2euler, euler2axangle, euler2mat
+from transforms3d.quaternions import quat2axangle, axangle2quat, mat2quat, qmult
+
+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 get_list_from_orix(grid,rounding=2):
     """
-    Creates a rotation list for the rotations within the fundamental zone of a given space group.
+    Converts an orix sample to a rotation list
 
     Parameters
     ----------
-    space_group_number : int
-        Between 1 and 230
-
-    resolution : float
-        The 'resolution' of the grid (degrees)
+    grid : orix.quaternion.rotation.Rotation
+        A grid of rotations
+    rounding : int, optional
+        The number of decimal places to retain, defaults to 2
 
     Returns
     -------
     rotation_list : list of tuples
+        A rotation list
     """
-    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)
+    z = grid.to_euler(convention="bunge")
+    rotation_list = z.data.tolist()
+    i = 0
+    while i < len(rotation_list):
+        rotation_list[i] = tuple(np.round(rotation_list[i],decimals=rounding))
+        i += 1
 
+    return rotation_list
 
-def get_grid_stereographic(crystal_system, resolution, equal="angle"):
+def get_fundamental_zone_grid(resolution=2, point_group=None,space_group=None):
     """
-    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.
+    Generates an equispaced grid of rotations within a fundamental zone.
 
     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.
+    resolution : float, optional
+        The characteristic distance between a rotation and its neighbour (degrees)
+    point_group : orix.quaternion.symmetry.Symmetry, optional
+        One of the 11 proper point groups, defaults to None
+    space_group: int, optional
+        Between 1 and 231, defaults to None
 
     Returns
     -------
     rotation_list : list of tuples
-        List of rotations
+        Grid of rotations lying within the specified fundamental zone
     """
-    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
 
+    orix_grid = get_sample_fundamental(resolution=resolution,space_group=space_group)
+    rotation_list = get_list_from_orix(orix_grid,rounding=2)
+    return rotation_list
 
-def get_local_grid(center, max_rotation, resolution):
+def get_local_grid(resolution=2, center=None, grid_width=10):
     """
-    Creates a rotation list for the rotations within max_rotation of center at a given rotation.
+    Generates a grid of rotations about a given rotation
 
     Parameters
     ----------
-    center : tuple
-        The orientation that acts as the center of the grid, as euler angles specified in the
-        'rzxz' convention (degrees)
+    resolution : float, optional
+        The characteristic distance between a rotation and its neighbour (degrees)
+    center : orix.quaternion.rotation.Rotation, optional
+        The rotation at which the grid is centered. If None (default) uses the identity
+    grid_width : float, optional
+        The largest angle of rotation away from center that is acceptable (degrees)
 
-    max_rotation : float
-        The largest rotation away from 'center' that should be included in the grid (degrees)
-
-    resolution : float
-        The 'resolution' of the grid (degrees)
+    See Also
+    --------
 
     Returns
     -------
     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)
-
-    return _returnable_eulers_from_axangle(raw_grid_axangle, "rzxz", round_to=2)
-
+    orix_grid =  get_sample_local(resolution=resolution, center=center, grid_width=grid_width)
+    rotation_list = get_list_from_orix(orix_grid,rounding=2)
+    return rotation_list
 
 def get_grid_around_beam_direction(beam_rotation, resolution, angular_range=(0, 360)):
     """
@@ -182,20 +146,91 @@ def get_grid_around_beam_direction(beam_rotation, resolution, angular_range=(0,
     axangle = euler2axangle(
         beam_rotation[0], beam_rotation[1], beam_rotation[2], "rzxz"
     )
-    euler_szxz = axangle2euler(axangle[0], axangle[1], "szxz")  # convert to szxz
-    rotation_alpha, rotation_beta = np.rad2deg(euler_szxz[0]), np.rad2deg(euler_szxz[1])
-
-    # see _create_advanced_linearly_spaced_array_in_rzxz for details
-    steps_gamma = int(np.ceil((angular_range[1] - angular_range[0]) / resolution))
-    alpha = np.asarray([rotation_alpha])
-    beta = np.asarray([rotation_beta])
-    gamma = np.linspace(
-        angular_range[0], angular_range[1], num=steps_gamma, endpoint=False
-    )
-    z = np.asarray(list(product(alpha, beta, gamma)))
-    raw_grid = Euler(
-        z, axis_convention="szxz"
-    )  # we make use of an uncommon euler angle set here for speed
-    grid_rzxz = raw_grid.to_AxAngle().to_Euler(axis_convention="rzxz")
-    rotation_list = grid_rzxz.to_rotation_list(round_to=2)
+    rotation_list = None
+    raise NotImplementedError("This functionality will be (re)added in future")
     return rotation_list
+
+
+def get_beam_directions_grid(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.
+    """
+    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

@@ -17,7 +17,6 @@
 # along with diffsims.  If not, see <http://www.gnu.org/licenses/>.
 
 import diffsims as ds
-from diffsims.generators.rotation_list_generators import get_grid_stereographic
 
 
 class StructureLibrary:
@@ -100,21 +99,20 @@ def from_crystal_systems(
             resolution in degrees
         equal : str
             Default is 'angle'
-        Returns
-        -------
-        StructureLibrary
+
+        Raises
+        ------
+        NotImplementedError:
+            "This function has been removed in version 0.3.0, in favour of creation from orientation lists"
         """
-        orientations = []
-        for system in systems:
-            orientations.append(get_grid_stereographic(system, resolution, equal))
-        return cls(identifiers, structures, orientations)
+        raise NotImplementedError("This function has been removed in version 0.3.0, in favour of creation from orientation lists")
 
     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/conftest.py

@@ -39,27 +39,3 @@ def default_simulator():
     accelerating_voltage = 300
     max_excitation_error = 1e-2
     return DiffractionGenerator(accelerating_voltage, max_excitation_error)
-
-
-@pytest.fixture()
-def random_eulers():
-    """ Using [0,360] [0,180] and [0,360] as ranges """
-    alpha = np.random.rand(100) * 360
-    beta = np.random.rand(100) * 180
-    gamma = np.random.rand(100) * 360
-    eulers = np.asarray((alpha, beta, gamma)).T
-    return eulers
-
-
-@pytest.fixture()
-def random_quats():
-    """ Unnormalised"""
-    q_rand = np.random.random(size=(1000, 4)) * 7
-    return q_rand
-
-
-@pytest.fixture()
-def random_axangles():
-    """ Unnormalised axes, & rotation between -pi and 2 pi """
-    axangle_rand = (np.random.random(size=(1000, 4)) * 3 * np.pi) - np.pi
-    return axangle_rand