diff --git a/diffsims/generators/rotation_list_generators.py b/diffsims/generators/rotation_list_generators.py index ea950be5..46d575c9 100644 --- a/diffsims/generators/rotation_list_generators.py +++ b/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,6 +60,47 @@ 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): """ @@ -65,7 +108,7 @@ def get_local_grid(center, max_rotation, resolution): 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") diff --git a/diffsims/tests/test_generators/test_rotation_list_generator.py b/diffsims/tests/test_generators/test_rotation_list_generator.py index e14b57a7..47ff027f 100644 --- a/diffsims/tests/test_generators/test_rotation_list_generator.py +++ b/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 diff --git a/diffsims/tests/test_utils/test_gridding_utils.py b/diffsims/tests/test_utils/test_gridding_utils.py index 450485ee..46d6c7fd 100644 --- a/diffsims/tests/test_utils/test_gridding_utils.py +++ b/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 diff --git a/diffsims/utils/gridding_utils.py b/diffsims/utils/gridding_utils.py index 47b71904..263e1882 100644 --- a/diffsims/utils/gridding_utils.py +++ b/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 diff --git a/diffsims/utils/vector_utils.py b/diffsims/utils/vector_utils.py index 60f1ceba..2264e8dd 100644 --- a/diffsims/utils/vector_utils.py +++ b/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