Add Sampling

.zenodo.json

@@ -34,6 +34,11 @@
       "name": "Anders Christian Mathisen",
       "affiliation": "Norwegian University of Science and Technology"
     },
+    {
+      "name": "Carter Francis",
+      "orcid": "0000-0003-2564-1851",
+      "affiliation": "University of Wisconsin Madison"
+    },
     {
       "name": "Simon Høgås"
     },

examples/rotations/sampling_rotations.py

@@ -0,0 +1,47 @@
+"""
+==================
+Sampling rotations
+==================
+
+This example shows how to sample some phase object in Orix. We will
+get both the zone axis and the reduced fundamental zone rotations for
+the phase of interest.
+"""
+
+from diffpy.structure import Atom, Lattice, Structure
+
+from orix.crystal_map import Phase
+from orix.sampling import get_sample_reduced_fundamental, get_sample_zone_axis
+from orix.vector import Vector3d
+
+a = 5.431
+latt = Lattice(a, a, a, 90, 90, 90)
+atom_list = []
+for coords in [[0, 0, 0], [0.5, 0, 0.5], [0, 0.5, 0.5], [0.5, 0.5, 0]]:
+    x, y, z = coords[0], coords[1], coords[2]
+    atom_list.append(Atom(atype="Si", xyz=[x, y, z], lattice=latt))  # Motif part A
+    atom_list.append(
+        Atom(atype="Si", xyz=[x + 0.25, y + 0.25, z + 0.25], lattice=latt)
+    )  # Motif part B
+struct = Structure(atoms=atom_list, lattice=latt)
+p = Phase(structure=struct, space_group=227)
+reduced_fun = get_sample_reduced_fundamental(resolution=4, point_group=p.point_group)
+
+vect_rot = (
+    reduced_fun * Vector3d.zvector()
+)  # get the vector representation of the rotations
+vect_rot.scatter(grid=True)  # plot the stereographic projection of the rotations
+
+# %%
+
+zone_axis_rot, directions = get_sample_zone_axis(
+    phase=p, density="7", return_directions=True
+)  # get the zone axis rotations
+zone_vect_rot = (
+    zone_axis_rot * Vector3d.zvector()
+)  # get the vector representation of the rotations
+zone_vect_rot.scatter(
+    grid=True, vector_labels=directions, text_kwargs={"size": 8, "rotation": 0}
+)  # plot the stereographic projection of the rotations
+
+# %%

orix/sampling/__init__.py

@@ -29,11 +29,18 @@
 )
 from orix.sampling.S2_sampling import sampling_methods as sample_S2_methods
 from orix.sampling.SO3_sampling import uniform_SO3_sample
-from orix.sampling.sample_generators import get_sample_fundamental, get_sample_local
+from orix.sampling.sample_generators import (
+    get_sample_fundamental,
+    get_sample_local,
+    get_sample_reduced_fundamental,
+    get_sample_zone_axis,
+)
 
 __all__ = [
     "get_sample_fundamental",
+    "get_sample_reduced_fundamental",
     "get_sample_local",
+    "get_sample_zone_axis",
     "sample_S2",
     "sample_S2_methods",
     "uniform_SO3_sample",

orix/sampling/sample_generators.py

@@ -20,10 +20,13 @@
 
 import numpy as np
 
-from orix.quaternion import OrientationRegion, Rotation, Symmetry
+from orix.crystal_map import Phase
+from orix.quaternion import OrientationRegion, Rotation, Symmetry, symmetry
 from orix.quaternion.symmetry import get_point_group
+from orix.sampling import sample_S2
 from orix.sampling.SO3_sampling import _three_uniform_samples_method, uniform_SO3_sample
 from orix.sampling._cubochoric_sampling import cubochoric_sampling
+from orix.vector import Vector3d
 
 
 def get_sample_fundamental(
@@ -165,3 +168,125 @@ def _remove_larger_than_angle(rot: Rotation, max_angle: Union[int, float]) -> Ro
     mask = half_angles < half_angle
     rot_out = rot[mask]
     return rot_out
+
+
+def get_sample_reduced_fundamental(
+    resolution: float,
+    mesh: str = None,
+    point_group: Symmetry = None,
+) -> Rotation:
+    """Produces rotations to align various crystallographic directions with
+    the z-axis, with the constraint that the first Euler angle phi_1=0.
+    The crystallographic directions sample the fundamental zone, representing
+    the smallest region of symmetrically unique directions of the relevant
+    crystal system or point group.
+    Parameters
+    ----------
+    resolution
+        An angle in degrees representing the maximum angular distance to a
+        first nearest neighbor grid point.
+    mesh
+        Type of meshing of the sphere that defines how the grid is created. See
+        orix.sampling.sample_S2 for all the options. A suitable default is
+        chosen depending on the crystal system.
+    point_group
+        Symmetry operations that determines the unique directions. Defaults to
+        no symmetry, which means sampling all 3D unit vectors.
+    Returns
+    -------
+    Rotation
+        (N, 3) array representing Euler angles for the different orientations
+    """
+    if point_group is None:
+        point_group = symmetry.C1
+
+    if mesh is None:
+        s2_auto_sampling_map = {
+            "triclinic": "icosahedral",
+            "monoclinic": "icosahedral",
+            "orthorhombic": "spherified_cube_edge",
+            "tetragonal": "spherified_cube_edge",
+            "cubic": "spherified_cube_edge",
+            "trigonal": "hexagonal",
+            "hexagonal": "hexagonal",
+        }
+        mesh = s2_auto_sampling_map[point_group.system]
+
+    s2_sample = sample_S2(resolution, method=mesh)
+    fundamental = s2_sample[s2_sample <= point_group.fundamental_sector]
+
+    phi = fundamental.polar
+    phi2 = (np.pi / 2 - fundamental.azimuth) % (2 * np.pi)
+    phi1 = np.zeros(phi2.shape[0])
+    euler_angles = np.vstack([phi1, phi, phi2]).T
+
+    return Rotation.from_euler(euler_angles, degrees=False)
+
+
+def _corners_to_centroid_and_edge_centers(corners):
+    """
+    Produces the midpoints and center of a trio of corners
+    Parameters
+    ----------
+    corners : list of lists
+        Three corners of a streographic triangle
+    Returns
+    -------
+    list_of_corners : list
+        Length 7, elements ca, cb, cc, mean, cab, cbc, cac where naming is such that
+        ca is the first corner of the input, and cab is the midpoint between
+        corner a and corner b.
+    """
+    ca, cb, cc = corners[0], corners[1], corners[2]
+    mean = tuple(np.add(np.add(ca, cb), cc))
+    cab = tuple(np.add(ca, cb))
+    cbc = tuple(np.add(cb, cc))
+    cac = tuple(np.add(ca, cc))
+    return [ca, cb, cc, mean, cab, cbc, cac]
+
+
+def get_sample_zone_axis(
+    density: str = "3",
+    phase: Phase = None,
+    return_directions: bool = False,
+) -> Rotation:
+    """Produces rotations to align various crystallographic directions with
+    the sample zone axes.
+
+    Parameters
+    ----------
+    density
+        Either '3' or '7' for the number of directions to return.
+    phase
+        The phase for which the zone axis rotations are required.
+    return_directions
+        If True, returns the directions as well as the rotations.
+    """
+    system = phase.point_group.system
+    corners_dict = {
+        "cubic": [(0, 0, 1), (1, 0, 1), (1, 1, 1)],
+        "hexagonal": [(0, 0, 1), (2, 1, 0), (1, 1, 0)],
+        "orthorhombic": [(0, 0, 1), (1, 0, 0), (0, 1, 0)],
+        "tetragonal": [(0, 0, 1), (1, 0, 0), (1, 1, 0)],
+        "trigonal": [(0, 0, 1), (-1, -2, 0), (1, -1, 0)],
+        "monoclinic": [(0, 0, 1), (0, 1, 0), (0, -1, 0)],
+    }
+    if density == "3":
+        direction_list = corners_dict[system]
+    elif density == "7":
+        direction_list = _corners_to_centroid_and_edge_centers(corners_dict[system])
+    else:
+        raise ValueError("Density must be either 3 or 7")
+
+    # rotate the directions to the z axis
+    rots = np.stack(
+        [
+            Rotation.from_align_vectors(v, Vector3d.zvector()).data
+            for v in direction_list
+        ]
+    )
+    rotations = Rotation(rots)
+    if return_directions:
+        return rotations, direction_list
+    else:
+        return rotations

orix/tests/sampling/test_sampling.py

@@ -16,12 +16,20 @@
 # You should have received a copy of the GNU General Public License
 # along with orix.  If not, see <http://www.gnu.org/licenses/>.
 
+from diffpy.structure import Atom, Lattice, Structure
 import numpy as np
 import pytest
 
+from orix.crystal_map import Phase
 from orix.quaternion import Rotation
-from orix.quaternion.symmetry import C2, C6, D6, get_point_group
-from orix.sampling import get_sample_fundamental, get_sample_local, uniform_SO3_sample
+from orix.quaternion.symmetry import C1, C2, C4, C6, D6, get_point_group
+from orix.sampling import (
+    get_sample_fundamental,
+    get_sample_local,
+    get_sample_reduced_fundamental,
+    get_sample_zone_axis,
+    uniform_SO3_sample,
+)
 from orix.sampling.SO3_sampling import _resolution_to_num_steps
 from orix.sampling._polyhedral_sampling import (
     _get_angles_between_nn_gridpoints,
@@ -175,3 +183,47 @@ def test_get_sample_fundamental_space_group(self, C6_sample):
         C2_sample = get_sample_fundamental(4, space_group=3, method="haar_euler")
         ratio = C2_sample.size / C6_sample.size
         assert np.isclose(ratio, 3, atol=0.2)
+
+    def test_get_sample_reduced_fundamental(self):
+        rotations = get_sample_reduced_fundamental(resolution=4)
+        rotations2 = get_sample_reduced_fundamental(resolution=4, point_group=C2)
+        rotations4 = get_sample_reduced_fundamental(resolution=4, point_group=C4)
+        rotations6 = get_sample_reduced_fundamental(resolution=4, point_group=C4)
+
+        assert (
+            np.abs(rotations.size / rotations2.size) - 2 < 0.1
+        )  # about 2 times more rotations
+        assert (
+            np.abs(rotations2.size / rotations4.size) - 2 < 0.1
+        )  # about 2 times more rotations
+        assert (
+            np.abs(rotations.size / rotations6.size) - 6 < 0.1
+        )  # about 6 times more rotations
+
+    @pytest.mark.parametrize("density", ("3", "7", "5"))
+    @pytest.mark.parametrize("get_directions", (True, False))
+    def test_get_zone_axis(self, density, get_directions):
+        a = 5.431
+        latt = Lattice(a, a, a, 90, 90, 90)
+        atom_list = []
+        for coords in [[0, 0, 0], [0.5, 0, 0.5], [0, 0.5, 0.5], [0.5, 0.5, 0]]:
+            x, y, z = coords[0], coords[1], coords[2]
+            atom_list.append(
+                Atom(atype="Si", xyz=[x, y, z], lattice=latt)
+            )  # Motif part A
+            atom_list.append(
+                Atom(atype="Si", xyz=[x + 0.25, y + 0.25, z + 0.25], lattice=latt)
+            )  # Motif part B
+        struct = Structure(atoms=atom_list, lattice=latt)
+        p = Phase(structure=struct, space_group=227)
+        if density == "5":
+            with pytest.raises(ValueError):
+                get_sample_zone_axis(phase=p, density=density)
+        else:
+            if get_directions:
+                rot, _ = get_sample_zone_axis(
+                    phase=p, density=density, return_directions=True
+                )
+            else:
+                rot = get_sample_zone_axis(phase=p, density=density)
+            assert isinstance(rot, Rotation)