improved testing, adressed comments PR

diffsims/generators/rotation_list_generators.py

@@ -22,7 +22,6 @@
 
 import numpy as np
 from scipy.spatial import cKDTree
-import math
 from itertools import product
 
 from orix.sampling.sample_generators import get_sample_fundamental, get_sample_local
@@ -161,7 +160,7 @@ def get_grid_around_beam_direction(beam_rotation, resolution, angular_range=(0,
     return rotation_list
 
 
-def normalize_vectors(vectors):
+def _normalize_vectors(vectors):
     """
     Helper function which returns a list of vectors normalized to length 1 from
     a 2D array representing a list of 3D vectors
@@ -171,7 +170,7 @@ def normalize_vectors(vectors):
 
 def get_uv_sphere_mesh_vertices(resolution):
     """
-    Return the vertices of a UV (polar) mesh on a unit sphere.
+    Return the vertices of a UV (spherical coordinate) mesh on a unit sphere.
     The mesh vertices are defined by the parametrization:
 
     x = sin(u)cos(v)
@@ -236,60 +235,72 @@ def get_cube_mesh_vertices(resolution, grid_type="spherified_corner"):
     Parameters
     ----------
     resolution : float
-        The maximum angle in degrees between neighboring grid points. 
-        Where this maximum angle lies depends on the grid_type.
-        To get an integer number of points between the cube face center
-        and the edges, the number of points is rounded up so that
-        resolution is always an upper limit.
+        The maximum angle in degrees between first nearest neighbor grid
+        points. 
     grid_type : str
         The type of cube grid, can be either `normalized` or `spherified_edge`
-        or `spherified_corner` (default).
-
-        In the normalized case, the grid on the surface of the cube is linear.
-        The maximum angle between nearest neighbors is found between the <001>
-        directions and the first grid point towards the <011> directions.
-
-        In the spherified_edge case, the grid is constructed so that there
-        are still two sets of perpendicular grid lines parallel to the {100}
-        directions on each cube face, but the spacing of the grid lines is
-        chosen so that the angles between the grid points on the line 
-        connecting the face centers (<001>) to the edges (<011>) are equal.
-        The maximum angle is also between the <001> directions and the first
-        grid point towards the <011> edges.
-
-        The spherified_corner case is similar to the spherified_edge case, but
-        the spacing of the grid lines is chosen so that the angles between
-        the grid points on the line connecting the face centers to the cube
-        corners (<111>) is equal. The maximum angle in this grid is from the
-        corners to the first grid point towards the cube face centers.
+        or `spherified_corner` (default). For details see notes.
 
     Returns
     -------
     points_in_cartesian : np.array (N,3)
         Rows are x, y, z where z is the 001 pole direction
 
+    Notes
+    -----
+    The resolution determines the maximum angle between first nearest
+    neighbor grid points, but to get an integer number of points between the
+    cube face center and the edges, the number of grid points is rounded up.
+    In practice this means that resolution is always an upper limit.
+    Additionally, where on the grid this maximum angle will be will depend
+    on the type of grid chosen. Resolution says something about the maximum
+    angle but nothing about the distribution of nearest neighbor angles or
+    the minimum angle - also this is fixed by the chosen grid.
+
+    In the normalized grid, the grid on the surface of the cube is linear.
+    The maximum angle between nearest neighbors is found between the <001>
+    directions and the first grid point towards the <011> directions. Points
+    approaching the edges and corners of the cube will have a smaller angular
+    deviation, so orientation space will be oversampled there compared to the
+    cube faces <001>.
+
+    In the spherified_edge grid, the grid is constructed so that there
+    are still two sets of perpendicular grid lines parallel to the {100}
+    directions on each cube face, but the spacing of the grid lines is
+    chosen so that the angles between the grid points on the line 
+    connecting the face centers (<001>) to the edges (<011>) are equal.
+    The maximum angle is also between the <001> directions and the first
+    grid point towards the <011> edges. This grid slightly oversamples the
+    directions between <011> and <111>
+
+    The spherified_corner case is similar to the spherified_edge case, but
+    the spacing of the grid lines is chosen so that the angles between
+    the grid points on the line connecting the face centers to the cube
+    corners (<111>) is equal. The maximum angle in this grid is from the
+    corners to the first grid point towards the cube face centers.
+
     References
     ----------
     [1] https://medium.com/game-dev-daily/four-ways-to-create-a-mesh-for-a-sphere-d7956b825db4
     """
     # the angle between 001 and 011
-    max_angle = np.arccos(1/np.sqrt(2))
+    max_angle = np.deg2rad(45)
     # the distance on the cube face from 001 to 011 grid point
     max_dist = 1
     if grid_type == "normalized":
         grid_len = np.tan(np.deg2rad(resolution))
-        steps = math.ceil(max_dist/grid_len)
+        steps = np.ceil(max_dist/grid_len)
         i = np.arange(-steps, steps)/steps
     elif grid_type == "spherified_edge":
-        steps = math.ceil(np.rad2deg(max_angle)/resolution)
+        steps = np.ceil(np.rad2deg(max_angle)/resolution)
         k = np.arange(-steps, steps)
         theta = np.arctan(max_dist)/steps
         i = np.tan(k*theta)
     elif grid_type == "spherified_corner":
         # the angle from 001 to 111
         max_angle_111 = np.arccos(1/np.sqrt(3))
         res_111 = np.deg2rad(resolution)
-        steps = math.ceil(max_angle_111/res_111)
+        steps = np.ceil(max_angle_111/res_111)
         k = np.arange(-steps, steps)
         theta = np.arctan(np.sqrt(2))/steps
         i = np.tan(k*theta)/np.sqrt(2)
@@ -312,13 +323,13 @@ def get_cube_mesh_vertices(resolution, grid_type="spherified_corner"):
                     [1, -1, -1]])
     # combine
     all_vecs = np.vstack([bottom, top, east, west, south, north, m_c])
-    return normalize_vectors(all_vecs)
+    return _normalize_vectors(all_vecs)
 
 
 def _compose_from_faces(corners, faces, n):
     """
     Helper function to refine a grid starting from a platonic solid,
-    adapted from meshzoo [1]
+    adapted from meshzoo
 
     Parameters
     ----------
@@ -334,6 +345,10 @@ def _compose_from_faces(corners, faces, n):
     -------
     vertices: numpy.array (N, 3)
         The coordinates of the refined mesh vertices.
+
+    See also
+    --------
+    :func:`get_icosahedral_mesh_vertices`
     """
     # create corner nodes
     vertices = [corners]
@@ -473,7 +488,7 @@ def _get_angles_between_nn_gridpoints(vertices, leaf_size=50):
     points on a grid of a sphere.
     """
     # normalize the vertex vectors
-    vertices = normalize_vectors(vertices)
+    vertices = _normalize_vectors(vertices)
     nn1_vec = _get_first_nearest_neighbors(vertices, leaf_size)
     # the dot product between each point and its nearest neighbor
     nn_dot = np.sum(vertices*nn1_vec, axis=1)
@@ -574,10 +589,30 @@ def get_icosahedral_mesh_vertices(resolution):
     return vertices
 
 
-def _vector_grid_to_euler(vectors):
+def beam_directions_grid_to_euler(vectors):
     """
-    Helper function to convert vectors representing zones to orientations
-    using only two Euler angles: (0, Phi, phi2).
+    Convert list of vectors representing zones to a list of Euler angles
+    in the bunge convention with the constraint that phi1=0.
+
+    Parameters
+    ----------
+    vectors: numpy.array (N, 3)
+        N 3-dimensional vectors to convert to Euler angles
+
+    Returns
+    -------
+    grid: numpy.array (N, 3)
+        Euler angles in bunge convention corresponding to each vector in
+        degrees.
+        
+    Notes
+    -----
+    The Euler angles represent the orientation of the crystal if that
+    particular vector were parallel to the beam direction [001]. The
+    additional constraint of phi1=0 means that this orientation is uniquely
+    defined for most vectors. phi1 represents the rotation of the crystal
+    around the beam direction and can be interpreted as the rotation of
+    a particular diffraction pattern.
     """
     norm = np.linalg.norm(vectors, axis=1)
     z_comp = vectors[:, 2]
@@ -617,7 +652,7 @@ def get_beam_directions_grid(crystal_system, resolution,
 
     mesh : str
         Type of meshing of the sphere that defines how the grid is created.
-        Options are: uv_sphere, normalized_cube, spherified_cube_corner 
+        Options are: uv_sphere, normalized_cube, spherified_cube_corner
         (default), spherified_cube_edge, and icosahedral.
 
     Returns
@@ -627,26 +662,21 @@ def get_beam_directions_grid(crystal_system, resolution,
     """
     if mesh == "uv_sphere":
         points_in_cartesians = get_uv_sphere_mesh_vertices(resolution)
-    elif mesh == "normalized_cube":
+    elif mesh == "normalized_cube" or "spherified_cube_edge":
         # special case: hexagon is a very small slice and 001 point can
         # be isolated. Hence we increase resolution to ensure minimum angle.
         if crystal_system == "hexagonal":
             resolution = np.rad2deg(
                     np.arctan(np.tan(np.deg2rad(resolution))/np.sqrt(2)))
-        points_in_cartesians = get_cube_mesh_vertices(resolution,
-                                                      grid_type="normalized")
-
-    elif mesh == "spherified_cube_edge":
-        # special case: hexagon is a very small slice and 001 point can
-        # be isolated.  Hence we increase resolution to ensure minimum angle.
-        if crystal_system == "hexagonal":
-            resolution = np.rad2deg(
-                    np.arctan(np.tan(np.deg2rad(resolution))/np.sqrt(2)))
-        points_in_cartesians = get_cube_mesh_vertices(resolution,
-                                                      grid_type="spherified_edge")
+        if mesh == "normalized_cube":
+            points_in_cartesians = get_cube_mesh_vertices(
+                    resolution, grid_type="normalized")
+        else:
+            points_in_cartesians = get_cube_mesh_vertices(
+                    resolution, grid_type="spherified_edge")
     elif mesh == "spherified_cube_corner":
-        points_in_cartesians = get_cube_mesh_vertices(resolution,
-                                                      grid_type="spherified_corner")
+        points_in_cartesians = get_cube_mesh_vertices(
+                resolution, grid_type="spherified_corner")
     elif mesh == "icosahedral":
         points_in_cartesians = get_icosahedral_mesh_vertices(resolution)
     else:
@@ -664,18 +694,18 @@ def get_beam_directions_grid(crystal_system, resolution,
         a, b, c = uvtw_to_uvw(a), uvtw_to_uvw(b), uvtw_to_uvw(c)
 
     # eliminates those points that lie outside of the stereographic triangle
-    sigfig = 7
+    epsilon = -1e13
     points_in_cartesians = points_in_cartesians[
-        np.round(np.dot(np.cross(a, b), c) * np.dot(np.cross(a, b),
-                 points_in_cartesians.T), sigfig) >= 0
+        np.dot(np.cross(a, b), c) *
+        np.dot(np.cross(a, b), points_in_cartesians.T) >= epsilon
     ]
     points_in_cartesians = points_in_cartesians[
-        np.round(np.dot(np.cross(b, c), a) * np.dot(np.cross(b, c),
-                 points_in_cartesians.T), sigfig) >= 0
+        np.dot(np.cross(b, c), a) *
+        np.dot(np.cross(b, c), points_in_cartesians.T) >= epsilon
     ]
     points_in_cartesians = points_in_cartesians[
-        np.round(np.dot(np.cross(c, a), b) * np.dot(np.cross(c, a),
-                 points_in_cartesians.T), sigfig) >= 0
+        np.dot(np.cross(c, a), b) *
+        np.dot(np.cross(c, a), points_in_cartesians.T) >= epsilon
     ]
 
     return points_in_cartesians

diffsims/tests/test_generators/test_rotation_list_generator.py

@@ -23,7 +23,14 @@
     get_grid_around_beam_direction,
     get_fundamental_zone_grid,
     get_beam_directions_grid,
-    _get_max_grid_angle
+    get_uv_sphere_mesh_vertices,
+    get_cube_mesh_vertices,
+    get_icosahedral_mesh_vertices,
+    beam_directions_grid_to_euler,
+    _normalize_vectors,
+    _get_first_nearest_neighbors,
+    _get_angles_between_nn_gridpoints,
+    _get_max_grid_angle,
 )
 
 
@@ -52,6 +59,88 @@ def test_get_grid_around_beam_direction():
     assert np.allclose([x[1] for x in grid], 90)  # taking z to y
 
 
+def test_get_uv_sphere_mesh_vertices():
+    grid = get_uv_sphere_mesh_vertices(10)
+    np.testing.assert_almost_equal(np.sum(grid), 0)
+    assert grid.shape[0] == 614
+    assert grid.shape[1] == 3
+    np.testing.assert_almost_equal(np.sum(grid), 0)
+    grid_unique = np.unique(grid, axis=0)
+    assert grid.shape[0] == grid_unique.shape[0]
+
+
+@pytest.mark.parametrize(
+    "grid_type,expected_len",
+    [
+        ("normalized", 866),
+        ("spherified_edge", 602),
+        ("spherified_corner", 866),
+    ]
+)
+def test_get_cube_mesh_vertices(grid_type, expected_len):
+    grid = get_cube_mesh_vertices(10, grid_type=grid_type)
+    assert grid.shape[0] == expected_len
+    assert grid.shape[1] == 3
+    np.testing.assert_almost_equal(np.sum(grid), 0)
+    grid_unique = np.unique(grid, axis=0)
+    assert grid.shape[0] == grid_unique.shape[0]
+    test_vectors = np.round(_normalize_vectors(np.array(
+                            [[1, 1, 1],
+                             [1, 0, 0],
+                             [1, 1, 0],
+                             [-1, 1, -1]])), 13).tolist()
+    grid = np.round(grid, 13)
+    for i in test_vectors:
+        assert i in grid.tolist()
+
+
+def test_get_cube_mesh_vertices_exception():
+    with pytest.raises(Exception):
+        get_cube_mesh_vertices(10, "non_existant")
+
+
+def test_first_nearest_neighbors():
+    grid = _normalize_vectors(np.array([[1, 0, 0],
+                                        [0, 1, 0],
+                                        [0, 1, 1],
+                                        [1, 0, 1],
+                                        ]))
+    fnn = _normalize_vectors(np.array([[1, 0, 1],
+                                       [0, 1, 1],
+                                       [0, 1, 0],
+                                       [1, 0, 0],
+                                       ]))
+    angles = np.array([45, 45, 45, 45])
+    fnn_test = _get_first_nearest_neighbors(grid)
+    angles_test = _get_angles_between_nn_gridpoints(grid)
+    assert np.allclose(fnn, fnn_test)
+    assert np.allclose(angles, angles_test)
+    assert _get_max_grid_angle(grid) == 45.
+
+
+def test_icosahedral_grid():
+    grid = get_icosahedral_mesh_vertices(10)
+    assert grid.shape[0] == 642
+    assert grid.shape[1] == 3
+    np.testing.assert_almost_equal(np.sum(grid), 0)
+    grid_unique = np.unique(grid, axis=0)
+    assert grid.shape[0] == grid_unique.shape[0]
+
+
+def test_vectors_to_euler():
+    grid = _normalize_vectors(np.array([[1, 0, 0],
+                                        [0, 1, 0],
+                                        [0, 1, 1],
+                                        [1, 0, 1],
+                                        ]))
+    ang = np.array([[0, 90, 0],
+                    [0, 90, 90],
+                    [0, 45, 90],
+                    [0, 45, 0],
+                    ])
+    assert np.allclose(ang, beam_directions_grid_to_euler(grid))
+
+
 @pytest.mark.parametrize(
     "mesh",
     [
@@ -78,3 +167,8 @@ def test_get_beam_directions_grid(crystal_system, mesh):
     grid = get_beam_directions_grid(crystal_system, 5, mesh=mesh)
     max_angle = _get_max_grid_angle(grid)
     assert max_angle <= 5
+
+
+def test_invalid_mesh_beam_directions():
+    with pytest.raises(Exception):
+        get_beam_directions_grid(10, "non_existant")