diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 647e4f60..1539477e 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -12,24 +12,34 @@ its best to adhere to `Semantic Versioning Added ----- +- ``Quaternion.mean(weight)`` allows the calculation of weighted averages. +- ``Misorientation.get_distance_matrix(lazy=True)``. + Currently opt-in, but will be the default in the next minor version. - Non-lazy computation of dot products with ``Misorientation.get_distance_matrix(lazy=True)``. Currently opt-in, but will be the default in the next minor version. Changed ------- +- Misorientations and Orientations account for symmetry when calculating means. - ``Miller.get_nearest()`` now raises a ``NotImplementedError`` rather than returning ``NotImplemented``. - ``Mille.mean(use_symmetry=True)`` now raises a ``NotImplementedError`` rather than returning ``NotImplemented``. - Improved (faster and using less memory) non-lazy computation of misorientation angles from ``Orientation.with_angle_outer()``. +- an OrientationRegion can now be calculated from symmetry for all 1024 possible combinations + of point groups (previously not implimented for 200 combinations). Fixed ----- - (Mis)orientation reduction to the fundamental zone via ``reduce()`` now correctly applies the symmetries in the opposite order, from right to left, `s_end * g * s_start`, where `g` is a (mis)orientation. +- Setting a Rotaion will now copy over the proper/improper marker if present. +- `Orientation.dot` and `Orientation.dot_outer` now correctly handle dot products for + multi-dimensional inputs. + 2026-06-06 - version 0.14.3 diff --git a/doc/user/bibliography.bib b/doc/user/bibliography.bib index 36ceb09a..c1259b28 100644 --- a/doc/user/bibliography.bib +++ b/doc/user/bibliography.bib @@ -77,6 +77,18 @@ @article{johnstone2020density volume = {53}, year = {2020} } +@article{markley_averaging_2007, + author = {Markley, F. Landis and Cheng, Yang and Crassidis, John L. and Oshman, Yaakov}, + doi = {10.2514/1.28949}, + journal = {Journal of Guidance, Control, and Dynamics}, + number = {4}, + title = {Averaging {Quaternions}}, + volume = {30}, + pages = {1193--1197}, + year = {2007}, + issn = {0731-5090, 1533-3884}, + url = {https://arc.aiaa.org/doi/10.2514/1.28949}, +} @PhdThesis{martineau2020multivariate, author = {Martineau, Benjamin Helks}, title = {Multivariate Analysis for Scanning (Transmission) Electron Diffraction}, diff --git a/examples/misorientations/reducing_and_averaging_misorientations.py b/examples/misorientations/reducing_and_averaging_misorientations.py new file mode 100644 index 00000000..7ee92131 --- /dev/null +++ b/examples/misorientations/reducing_and_averaging_misorientations.py @@ -0,0 +1,230 @@ +# Copyright 2018-2026 the orix developers +# +# This file is part of orix. +# +# orix is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# orix is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with orix. If not, see . +# +r""" +======================================================= +Reducing and Averaging Misorientations and Orientations +======================================================= + +This example introduces the concept of reducing an orientation or +misorientation with respect to symmetry, as well as the related +concept of averaging a misorientation. +""" + +import matplotlib.pyplot as plt +import numpy as np +from scipy.stats import norm + +from orix.plot.rotation_plot import _setup_rotation_plot +import orix.quaternion as oqu +import orix.quaternion.symmetry as osm + +# Reproducible random data. +np.random.seed(2319) + + +# convenience function for plotting. +def setup_plot(fz): + fig, ax = _setup_rotation_plot(projection="homochoric") + fig.set_figwidth(6) + fig.set_figheight(5) + ax.axis("off") + ax._correct_aspect_ratio(fz) + return (fig, ax) + + +######################################################################################### +# Basic Example +# ------------- +# `Orientation.reduce` and `Misorientation.reduce` will convert a transform +# with symmetry to a unique equivalent representation with the smallest +# possible angle of rotation. This is what is meant by "reducing" a +# transform. + +g_random = oqu.Orientation.random(shape=1, symmetry=osm.Oh) +g_reduced = g_random.reduce() + +g_fz = oqu.OrientationRegion.from_symmetry(osm.C1, osm.Oh) +g_all = osm.Oh.outer(g_random) +g_equiv = g_all[oqu.Rotation(g_all).dot(g_reduced) < 0.999] + +fig_1, ax_1 = setup_plot(g_fz) +ax_1.set_xlim([-1, 1]) +ax_1.set_aspect("equal") +ax_1.scatter(g_equiv, color="black") +ax_1.scatter(g_reduced, color="red") +ax_1.plot_wireframe(g_fz, color="grey") +fig_1.suptitle( + "Reduced (red) and equivalent (black) \nrepresentations in point group Oh (m-3m)" +) + +######################################################################################### +# The inclusion of symmetry combined with the periodic nature of rotations +# can make the definition of a mean or average ambiguous (more on this below), +# so the first step when an average is calculated is to reduce the transforms +# and then calculate the average. A side effect of this is the mean +# (mis)orientation returned by ORIX will also be a reduced representation. + +qu_data = np.stack([norm.rvs(i, 0.06, 20) for i in [0.39, 0.28, -0.39, 0.78]]).T +m_sym = [osm.C3, osm.D2] + +m_clustered = oqu.Misorientation(qu_data, symmetry=m_sym) +m_reduced = m_clustered.reduce() +m_mean, m_neighbors = m_clustered.mean(return_neighbors=True) + +fz_cluster = oqu.OrientationRegion.from_symmetry(*m_sym) + +fig2, ax2 = setup_plot(fz_cluster) +ax2.scatter(oqu.Rotation(m_clustered), c="k") +ax2.scatter(oqu.Rotation(m_reduced), c="r") +ax2.scatter(oqu.Rotation(m_mean), color="blue", marker="X", s=100) +ax2.plot_wireframe(fz_cluster, color="grey") +fig2.suptitle( + "Reduced (red), Original (black), and symmetry-aware \nMean (blue) for {C3(3)-> D2(222)} system" +) + +######################################################################################### +# The Fundamental Zone +# -------------------- +# +# In the plots above, wireframes were included that defined bounded volumes +# within which all reduced representations fell. This is known as a +# Fundamental Zone (FZ), and contains a single unique (aka, fundamental) +# representation of every transformation with respect to a given symmetry. +# representations that fall within a fundamental zone are also garunteed to +# have the smallest possible angular component. +# +# There are several ways in which a fundamental zone can be defined, with most +# discrepencies stemming from how improper transforms +# (ie, inversions and rotoinversions) should be handled. ORIX uses the rules +# presented in :cite:`krakow2017onthree`, but expanded to all 1024 +# misorientation groups. This can be verified by comparing the following plots +# to Figure 5 of the same paper. + +name2group = {x.name: x for x in osm._groups} + +fig5 = plt.figure(figsize=[5, 8]) + +base_pairs = [ + ["432", "3"], + ["23", "3"], + ["432", "1"], + ["23", "1"], + ["3", "422"], + ["622", "1"], + ["4", "211"], + ["32", "1"], + ["222", "1"], + ["3", "4"], + ["6", "1"], + ["4", "1"], + ["3", "1"], + ["211", "1"], + ["1", "1"], +] +for i, pair in enumerate(base_pairs): + s1 = name2group[pair[0]] + s2 = name2group[pair[1]] + fz = oqu.OrientationRegion.from_symmetry(s1, s2) + ax = fig5.add_subplot(5, 3, i + 1, projection="homochoric") + ax.axis("off") + ax._correct_aspect_ratio(fz) + ax.plot_wireframe(fz) + ax.set_title("{}: {} -> {}".format("abcdefghijklmno"[i], s1.name, s2.name)) + +plt.tight_layout() + +######################################################################################### +# If users are unfamiliar with how Rodrigues, Homochoric, or NeoEulerian plots +# are used to plot rotations in 3D space, the same paper also contains a +# concise overview. +# +# It is not ran as part of this example since calculating and plotting 66 +# fundamental zones is time consuming, but the following plot matches Table +# 3 of the same paper, showing the subdivisions of the above 15 zones caused +# by shared rotation elements between the point groups. Users may plot it +# themselves by running this code locally with `plot_me=True` + + +plot_me = False + +names = ["432", "23", "622", "6", "32", "3", "422", "4", "222", "211", "1"] +if plot_me is True: + table_fig = plt.figure() + for i in range(len(names)): + n1 = names[i] + s1 = name2group[n1] + for j in range(len(names) - i): + n2 = names[-j - 1] + s2 = name2group[n2] + fz = oqu.OrientationRegion.from_symmetry(s1, s2) + ij = j * 11 + i + 1 + ax = table_fig.add_subplot(11, 11, ij, projection="homochoric") + ax.axis("off") + ax._correct_aspect_ratio(fz) + ax.plot_wireframe(fz) + +######################################################################################### +# Defining a Mean in Rotation Space +# --------------------------------- +# +# Up until now, we have not defined what is meant by the mean of a group of +# transforms. Because rotation space is periodic, the concept of a Euclidean +# norm does not apply, and instead a Frobenius norm is used, which in this +# context can be thought of as the magnitude of the angular rotation necessary +# to align two transforms. Noteably, this is NOT simply the normalized average +# of two transform's quaternion representations, as is sometimes done in other +# software to get a fast approximation for clustered transforms. See the +# docstring for `Quaternion.mean` for details on this topic. +# +# The extension of this is a mean transform is defined as the transform whose +# total angular deviation from all transforms in the queried group is the +# minimum possible value. +# +# However, there is not a convenient algorithm to calculate this correctly, +# so instead ORIX does the following: +# +# 1) transforms are reduced to the appropriate fundamental zone. +# 2) A rough mean is calculated. +# 3) transforms with equivalents closer to the rough mean are updated to the closer value +# 4) A precise mean is recalculated. +# +# The plot below is provided to help visualize this process. + +np.random.seed(2319) +qu_data = np.stack([norm.rvs(i, 0.1, 20) for i in [0.1, 0.1, 0.2, 0.3]]).T + +o_cluster = oqu.Orientation(qu_data, symmetry=osm.D2) +o_reduced = o_cluster.reduce() +o_mean, o_neighbors = o_cluster.mean(return_neighbors=True) +o_flipped = o_neighbors[np.abs(o_neighbors.angle - o_reduced.angle) > 1e-3] + + +fig_ave, ax_ave = setup_plot(fz) +ax_ave.scatter(oqu.Rotation(o_cluster), color="black") +ax_ave.scatter(oqu.Rotation(o_reduced), color="red") +ax_ave.scatter(oqu.Rotation(o_flipped), color="green") +ax_ave.scatter(oqu.Rotation(o_mean), color="blue", marker="X", s=100) + +fz = oqu.OrientationRegion.from_symmetry(end=osm.D2) +ax_ave.plot_wireframe(fz, color="grey") +ax_ave.set_xlim([-1, 1]) +ax_ave.set_ylim([-1, 1]) +ax_ave.set_zlim([-1, 1]) +fig_ave.suptitle( + "Original (black), Reduced (red), and flipped(green) \nrepresentations for point group D2 (222)" +) diff --git a/examples/plotting/visualizing_rotation_vector_paths.py b/examples/plotting/visualizing_rotation_vector_paths.py index eb92c998..84f38318 100644 --- a/examples/plotting/visualizing_rotation_vector_paths.py +++ b/examples/plotting/visualizing_rotation_vector_paths.py @@ -103,11 +103,11 @@ path.symmetry = Oh colors2 = mpl.colormaps["inferno"](np.linspace(0, 1, n_steps)) -fig = path.scatter("rodrigues", position=121, return_figure=True, c=colors2) -path.scatter("ipf", position=122, figure=fig, c=colors2) +fig2 = path.scatter("rodrigues", position=121, return_figure=True, c=colors2) +path.scatter("ipf", position=122, figure=fig2, c=colors2) # Plot the rest -rod_ax, ipf_ax = fig.axes +rod_ax, ipf_ax = fig2.axes rod_ax.set_title("Orientation paths in Rodrigues space") ipf_ax.set_title("Vector paths in IPF-Z", pad=15) @@ -125,7 +125,9 @@ # # Rotate vectors around the (1, 1, 1) axis on a stereographic plot. -vec_ax = plt.subplot(projection="stereographic") + +fig3 = plt.figure() +vec_ax = fig3.add_subplot(111, projection="stereographic") vec_ax.set_title(r"Stereographic") vec_ax.set_labels("X", "Y") diff --git a/orix/quaternion/misorientation.py b/orix/quaternion/misorientation.py index fbaa8698..b085b99d 100644 --- a/orix/quaternion/misorientation.py +++ b/orix/quaternion/misorientation.py @@ -31,7 +31,6 @@ from scipy.spatial.transform import Rotation as SciPyRotation from tqdm import tqdm -from orix._utils.deprecation import deprecated from orix.quaternion._numba import _mori_distance_matrix from orix.quaternion.orientation_region import OrientationRegion from orix.quaternion.rotation import Rotation @@ -393,37 +392,13 @@ def equivalent(self, grain_exchange: bool = False) -> Misorientation: return self.__class__(equivalent).flatten() - @deprecated(since="0.14", removal="0.15", alternative="reduce") - def map_into_symmetry_reduced_zone(self, verbose: bool = False) -> Misorientation: - """Return equivalent transformations which have the smallest - angle of rotation as a new misorientation. - - Parameters - ---------- - verbose - Whether to print a progressbar. Default is ``False``. - - Returns - ------- - M - A new misorientation object with the assigned symmetry. - - Examples - -------- - >>> from orix.quaternion.symmetry import C4, C2 - >>> data = np.array([[0.5, 0.5, 0.5, 0.5], [0, 1, 0, 0]]) - >>> M = Misorientation(data) - >>> M.symmetry = (C4, C2) - >>> M.map_into_symmetry_reduced_zone() - Misorientation (2,) 4, 2 - [[-0.7071 0.7071 0. 0. ] - [ 0. 1. 0. 0. ]] - """ - return self.reduce(verbose) - def reduce(self, verbose: bool = False) -> Misorientation: - """Return equivalent transformations which have the smallest - angle of rotation as a new misorientation. + """Return symmetrically equivalent transformations with the + smallest angle of rotation. + + For misorientations, reduced representations are further + restricted to transforms inside their symmetry's fundamental + zone. See Notes section for details. Parameters ---------- @@ -448,15 +423,32 @@ def reduce(self, verbose: bool = False) -> Misorientation: Notes ----- - The misorientation with the smallest rotation angle is the one - inside the misorientation fundamental zone (FZ, asymmetric - domain) given by the proper point groups of :attr:`symmetry`. - The definition of the FZ follows the procedure in section 5.3.1 - in :cite:`martineau2020multivariate`. - - An alternative description of finding the Rodrigues FZ is given - in :cite:`morawiec1996rodrigues`, which is the basis for - reduction of (mis)orientations in EMsoft. + In ORIX, fundamental zones are defined as bounded volumes in + quaternion space containing only proper rotations. This is + a common convention, for reasons discussed in the docstring + of :func:`orix.quaternion.OrientationRegion.from_symmetry()`. + + This is relevant for defining a reduced representation + because the brute force expansion of a misorientation to it's + symmetric equivalents can produce up to four unique + rotations that fall inside the fundamental zone and also all + have the same rotation angle. + + For all orientations and 924 of the possible 1024 + misorientation symmetries, this fact is irrelevant, as either + only a single proper rotation or an identical proper + and improper pair will map to the fundamental zone. + + The remaining 100 cases occur when both symmetries are + improper and don't possess an inversion. In this case, four + unique values fall within the fundamental zone, including + two proper rotations. One of these is produced using only + proper rotations, whereas the second is a pseudo-proper + rotatoin resulting from two consecutive rotoinversions. + + Since pseudo-proper variants cannot be reached through any + combination of proper rotations from either symmetry, they + are ignored by ORIX and only the proper rotation is returned. """ # Combine symmetry elements of start and end of transformation # given by the (mis)orientation @@ -469,16 +461,21 @@ def reduce(self, verbose: bool = False) -> Misorientation: # (orientation) or MacKenzie (misorientation) fundamental zone # (FZ), given by the symmetry elements. We loop over all # symmetry pairs and rotate all (mis)orientations until all are - # inside the FZ. - fz = OrientationRegion.from_symmetry(s1=start, s2=end) + # inside the FZ. Ignore symmetry combinations that need an inversion, + # as these are not handled by the MacKenzie/Rodrigues definitions of + # a fundamental zone. + fz = OrientationRegion.from_symmetry(start=start, end=end) reduced = self.__class__.identity(self.shape) is_outside = np.ones(self.shape, dtype=bool) for sym_start, sym_end in symmetry_pairs: + if sym_start.improper != sym_end.improper: + continue reduced[is_outside] = sym_end * self[is_outside] * sym_start is_outside = ~(reduced < fz) if not is_outside.any(): break - + # convert to northern hemisphere representations + reduced.data[reduced.a < 0] = reduced.data[reduced.a < 0] * -1 reduced._symmetry = (start, end) return reduced @@ -672,6 +669,117 @@ def inv(self) -> Misorientation: r"""Return the inverse misorientations :math:`M^{-1}`.""" return self.__invert__() + def mean( + self, + weights: np.ndarray | None = None, + include_improper: bool = False, + ignore_symmetry: bool = False, + return_neighbors: bool = False, + verbose: bool = False, + ) -> Misorientation: + """Return the symmetry-respecting mean (mis)orientation. + + Parameters + ---------- + weights + An optional array of weights for calculating a weighted + average instead of the unweighted mean. Must be the same + size as the quaternion array. + + include_improper + If True, equivalent representations that require inversion + symmetry to calculate will be excluded. See Notes for + details. Default is False. + + ignore_symmetry + If True, ignore all symmetry considerations. See Notes + for detials. Default is False. + + return_neighbors + If True, returns the nearest neighbors used to calculate + the mean. Default is False. + + verbose + If True, print progress bars. Default is False. + + Returns + ------- + mean + Mean (mis)orientation. + + neighbors + If `return_neighbors` is True, returns the + representations used to calculate the mean. + + Notes + ----- + This method uses the Frobenius norm of rotation space to + define a mean for rotations, as given in equations 12 and 13 + of :cite:`markley_averaging_2007`. Refer to + :func:`orix.quaternion.Quaternion.mean` for details. + + To account for symmetry, the following proceedure is used: + + 1) Misorientations are reduced to the fundamental zone. + 2) The rough mean is calculated. + 3) Misorientations with equivalent values closer to the + rough mean are updated to the nearby value. + 4) The precise mean is recalculated. + + if ``ignore_symmetry`` is True, steps 3 and 4 are skipped, + and the mean is given as a Rotation to signify the loss of + symmetry information. + + Since a pure rotation cannot align an inverted reference + frame with an uninverted one, a Frobenius norm cannot be + calculated for a mix of proper and improper rotations. + By default, this problem is addressed by ignoring + symmetrically equivalent operations that include inversion. + This aligns with the definition of a fundamental zone in + orientation space used in + :func:`orix.quaternion.Misorientation.reduce` and + :func:`orix.quaternion.OrientationRegion.from_symmetry`. + Setting `include_improper=True` will instead investigate all + symmetry options and treat all options as proper rotations + when calculating the mean. + """ + if ignore_symmetry is True: + # convert to a rotation to emphasize loss of symmetry information + return Rotation(self.data).mean(weights=weights) + + if verbose: + print("reducing to fundamental zone...") + # overwrite new nearest values into neighbors. Use rot for calculating + # candidate for inclusion in neighbors. + neighbors = self.reduce(verbose=verbose) + rots = Rotation(neighbors.data) + rough_mean = rots.mean(weights=weights) + + max_dp = rots.dot(rough_mean) + start, end = self._symmetry + if not include_improper: + start = start.proper_subgroup + end = end.proper_subgroup + symmetry_pairs = iproduct(Rotation(start), Rotation(end)) + if verbose: + print("checking for closer equivalent representations...") + s = start.size * end.size + symmetry_pairs = tqdm(symmetry_pairs, total=s) + for start, end in symmetry_pairs: + candidates = end * rots * start + dp = np.abs(candidates.dot(rough_mean)) + mask = dp > max_dp + # copy quaternion plus improper marker + neighbors._data[mask, :] = candidates._data[mask, :] + max_dp[dp > max_dp] = dp[dp > max_dp] + + fine_mean_rot = Rotation(neighbors.data).mean(weights=weights) + fine_mean = self.__class__(fine_mean_rot) + fine_mean._symmetry = self._symmetry + if return_neighbors: + return [fine_mean, neighbors] + return fine_mean + def _get_distance_matrix_dask( mori: Misorientation, diff --git a/orix/quaternion/orientation.py b/orix/quaternion/orientation.py index 1e5319a6..5bdc59e9 100644 --- a/orix/quaternion/orientation.py +++ b/orix/quaternion/orientation.py @@ -638,11 +638,9 @@ def dot(self, other: Orientation) -> np.ndarray: >>> O1.dot(O2) array([0.92387953, 0.92387953]) """ - symmetry = _get_unique_symmetry_elements(self.symmetry, other.symmetry) - M = other * ~self - all_dot_products = Rotation(M).dot_outer(symmetry) - highest_dot_products = np.max(all_dot_products, axis=-1) - return highest_dot_products + sym = [self.symmetry, other.symmetry] + M = Misorientation(other * ~self, symmetry=sym).reduce() + return M.a def dot_outer(self, other: Orientation) -> np.ndarray: """Return the symmetry reduced dot products of all orientations @@ -671,15 +669,9 @@ def dot_outer(self, other: Orientation) -> np.ndarray: array([[0.92387953, 1. ], [1. , 0.92387953]]) """ - symmetry = _get_unique_symmetry_elements(self.symmetry, other.symmetry) - M = other.outer(~self) - all_dot_products = Rotation(M).dot_outer(symmetry) - highest_dot_products = np.max(all_dot_products, axis=-1) - # need to return axes order so that self is first - order = tuple(range(self.ndim, self.ndim + other.ndim)) + tuple( - range(self.ndim) - ) - return highest_dot_products.transpose(*order) + sym = [self.symmetry, other.symmetry] + M = Misorientation(other.outer(~self), symmetry=sym).reduce() + return M.a def plot_unit_cell( self, diff --git a/orix/quaternion/orientation_region.py b/orix/quaternion/orientation_region.py index 64ca72ca..317309fe 100644 --- a/orix/quaternion/orientation_region.py +++ b/orix/quaternion/orientation_region.py @@ -24,14 +24,42 @@ import numpy as np from orix._utils import constants +from orix._utils.deprecation import deprecated, deprecated_argument from orix.quaternion.quaternion import Quaternion from orix.quaternion.rotation import Rotation from orix.quaternion.symmetry import C1, Symmetry, get_distinguished_points from orix.vector.neo_euler import Rodrigues -def _get_large_cell_normals(s1, s2): - dp = get_distinguished_points(s1, s2) +def _get_large_cell_normals( + start: Symmetry = C1, end: Symmetry = C1, dp: Rotation = None +) -> Rotation: + """Return rotations defining fundamental zone bounds due to + symmetry. + + Given two symmetries, calculates every unique rotation equivalent + to the identity, called "distinguished points". A Voronoi + tessletation is then done to define the bounds within which all + rotations are closer to identity than any distinguished points. + + Instead of calculating them, a set of distinguished points can + also be added by the user, in which case the symmetry arguments + will be ignored. This is useful for defining non-crystallographic + orientation regions. + Parameters + ---------- + start + Starting symmetry + end + Ending symmetry + dp + Optional. set of rotations that are equivalent to identity. + If given, s1 and s2 will be ignored, and the normals will + be defined via tesselation of these points. + + """ + if dp is None: + dp = get_distinguished_points(start, end) if dp.size == 0: return Rotation.empty() @@ -56,71 +84,99 @@ def _get_large_cell_normals(s1, s2): return normals -def get_proper_groups(Gl: Symmetry, Gr: Symmetry) -> tuple[Symmetry, Symmetry]: +@deprecated(since="0.15", removal="0.16") +def get_proper_groups(start: Symmetry, end: Symmetry) -> tuple[Symmetry, Symmetry]: """Return the appropriate groups for the asymmetric domain calculation. Parameters ---------- - Gl - First point group. - Gr - Second point group. + start + Initial point group, C1 for orientations. + end + Ending point group. Returns ------- - Gl - First proper subgroup(s) or proper inversion subgroup(s), as - appropriate. - Gr - Second proper subgroup(s) or proper inversion subgroup(s), as - appropriate. - - Raises - ------ - NotImplementedError - If both groups are improper and neither contain an inversion, - special consideration is needed which is not yet implemented in - orix. + start + Initial proper, inversion, or laue subgroup as appropriate. + end + Final proper, inversion, or laue subgroup as appropriate. + + Notes + ----- + ORIX follows the asymmetric domain/fundamental zone definitions + from in :cite:`krakow2017onthree`. However, for reasons given + in section 3(b), that paper only defines domains for the 121 + combinations of proper point groups. Orix extends their logic for + defining domains to the remaining 903 point group combinations, + which are always repetitions of the original 121. + + This expansion is unimportant for all orientations as well + as any misorientations where at least on point group is proper + and/or centrosymmetric, which together accounts for 924 of the + possible 1024 possible misorientation symmetries. This includes + all data from EBSD due to the artificial centrosymmetry caused by + kikuchi diffraction. For the remaining 100 misorientations where + both point groups are improper but don't contain an inversion + (for example, 6mm-->6mm), there are up to four rotations that + map to the fundamental zone; a proper, two improper, and a + pseudo-proper created from two roto-inversions. Common practice + in these edge cases is to ignore all but the proper rotation, as + is done in :func:`Misorientation.reduce`. """ - if Gl.is_proper and Gr.is_proper: - return Gl, Gr - elif Gl.is_proper and not Gr.is_proper: - return Gl, Gr.proper_subgroup - elif not Gl.is_proper and Gr.is_proper: - return Gl.proper_subgroup, Gr + if start.is_proper and end.is_proper: + return start, end + elif start.is_proper and not end.is_proper: + return start, end.proper_subgroup + elif not start.is_proper and end.is_proper: + return start.proper_subgroup, end else: - if Gl.contains_inversion and Gr.contains_inversion: - return Gl.proper_subgroup, Gr.proper_subgroup - elif Gl.contains_inversion and not Gr.contains_inversion: - return Gl.proper_subgroup, Gr.laue_proper_subgroup - elif not Gl.contains_inversion and Gr.contains_inversion: - return Gl.laue_proper_subgroup, Gr.proper_subgroup + if start.contains_inversion and end.contains_inversion: + return start.proper_subgroup, end.proper_subgroup + elif start.contains_inversion and not end.contains_inversion: + return start.proper_subgroup, end.laue_proper_subgroup + elif not start.contains_inversion and end.contains_inversion: + return start.laue_proper_subgroup, end.proper_subgroup else: - raise NotImplementedError( - "Both groups are improper, " "and do not contain inversion." - ) + return start.laue_proper_subgroup, end.laue_proper_subgroup class OrientationRegion(Rotation): - """Some subset of the complete space of orientations. - - The complete orientation space represents every possible orientation - of an object. The whole space is not always needed, for example if - the orientation of an object is constrained or (most commonly) if - the object is symmetrical. In this case, the space can be segmented - using sets of Rotations representing boundaries in the space. This - is clearest in the Rodrigues parametrisation, where the boundaries - are planes, such as the example here: the asymmetric domain of an - adjusted 432 symmetry. + """A subset of rotation space. + + The complete set of all possible rigid body rotations is called + SO(3). it can be thought of as half the quaternion unit sphere, + the entirety of Rodrigues space, the set of all 3x3 matrices + with a determinate of 1, or various other descriptors based + on the application. + + Sometimes, this whole space is not need, for example if the + orientation of an object is constrained or (most commonly) if + the object is symmetrical. In this case, the space can be + segmented using set of rotations representing boundaries in the + space. This can be most easily visualized using Rodrigues + space, where the boundaries become flat planes normal to the + rodrigues vectors of those bounding rotations. .. image:: /_static/img/orientation-region-Oq.png :width: 300px :alt: Boundaries of an orientation region in Rodrigues space. :align: center - Rotations or orientations can be inside or outside of an orientation - region. + Quaternions can then be quickly defined as inside or outside of + these regions via a dot product operation. + + Notably, these regions are only defined in SO(3), which means + they cannot account for improper operations. This is why + OrientationRegion.from_symmetry() calculates identical regions + for point groups 432 and m-3m despite m-3m having twice as many + distinguished points. This ends up being irrelevant for + Orientations since any improper operations that place a point + within a fundamental zone always have a paired proper operation + that returns an identical quaternion, but it can create confusion + for misorientations with rotoinversions when users assume an + OrientationRegion can uniquely define a true fundamental zone. """ # ------------------------ Dunder methods ------------------------ # @@ -128,8 +184,9 @@ class OrientationRegion(Rotation): def __gt__(self, other: OrientationRegion) -> np.ndarray: """Overridden greater than method. - Applying this to an orientation will return only those that lie - within the region. + Applying this to a rotation will return only those that lie + within the region. This operation does not account for + inversion. """ c = Quaternion(self).dot_outer(Quaternion(other)) inside = np.logical_or( @@ -141,25 +198,115 @@ def __gt__(self, other: OrientationRegion) -> np.ndarray: # ------------------------ Class methods ------------------------- # @classmethod - def from_symmetry(cls, s1: Symmetry, s2: Symmetry = C1) -> OrientationRegion: - """Return the set of unique (mis)orientations of a symmetrical - object. + @deprecated_argument("s1", since="0.15", removal="0.16", alternative="start") + @deprecated_argument("s2", since="0.15", removal="0.16", alternative="end") + def from_symmetry( + cls, + start: Symmetry = C1, + end: Symmetry = C1, + s1=None, + s2=None, + ) -> OrientationRegion: + """Return an orientation region for a given symmetry. + + These regions are identical to the fundamental zone for all + orientations and every misorientation where both + symmetries are proper and/or centrosymmetric. For + all other cases, it is still garunteed to inlude only one + unique represenation achievable though proper rotations.See + Notes for details. Parameters ---------- - s1 - First symmetry. - s2 - Second symmetry. Default is C1 (the identity). + start + Initial point group, C1 for passive orientations. + end + Ending point group. Returns ------- region The orientation region. + + Notes + ----- + ORIX follows the asymmetric domain/fundamental zone definitions + from in :cite:`krakow2017onthree`. However, for reasons given + in section 3(b), domains are only described for the 121 + combinations of proper point groups. Section 5.3.1 of + :cite"martineau2020multivariate" gives pseudocode for + extending this original logic to define domains for all 1024 + possible symmetry cases. + + This ends up being a trivial terminology issue for all + orientations as well as any misorientations where both point + groups are proper and/or centrosymmetric. This includes all + EBSD data as well, since kikuchi diffraction introduces an + artificial centrosymmetry. For these 704 cases, either the + region returned bounds a fully unique zone, or it bounds + all proper representations, and the improper representations + have identical quaternion representations. + + For the remaining 320 misorientations where one or both + point groups contain rotoinversions but are not + centrosymmetric (for example, 6mm --> 6mm), there are always + one and possibly two improper rotations that also map to the + orientation region but with unique quaternion values, as well + as a possible unique pseudo-proper rotation only achievable + through two rotoinversions. There is currently no concensus + on how to define unique fundamental zones for these edge + cases. """ - s1, s2 = get_proper_groups(s1, s2) - large_cell_normals = _get_large_cell_normals(s1, s2) - disjoint = s1 & s2 + if s1 is not None: + start = s1 + if s2 is not None: + end = s2 + # Step 1: fundamental zones are only defined for proper rotations. + # add inversion centers where necessary to define as unique as + # possible of a fundamental zone, then remove all improper operators. + + # If either symmetry is proper, any improper symmetries of the second + # group will fall outside the shared fundamental zone. + if not start.is_proper and not end.is_proper: + # If both symmetries contain an inversion, all improper operators + # will have a paired proper operator. If neither do, the proper + # and improper rotations will form two identical but inverted + # fundamental zones. Both cases produce one proper, two improper, + # and one pseudo-proper fundamental zone, but the first case is + # uninteresting because they have idencial quaternion + # representations. The second case requires some special + # consideration when using fundamental zones for averages or + # reducing. Regardless though, for both cases, reducing both + # symmetries to their proper subgroup gives the correct + # fundamental zone. + if start.contains_inversion != end.contains_inversion: + # The remaining case is when only one of the two groups + # contains an inversion. In this case, the combination of + # an inversion and rotation creates a mirror that needs to + # be added to the disjoint group. this is easiest done + # by converting the inversion-less symmetry to it's laue + # symmetry. + if not start.contains_inversion: + start = start.laue + if not end.contains_inversion: + end = end.laue + # With mirrors from inversion/rotoinversion combinations accounted + # for, remove all improper operators. The fundamental zone will now be + # one of the 61 in Table 3 of krakow2017. + start = start.proper_subgroup + end = end.proper_subgroup + + # Step 2: define the bounding cells using the distinguished points. + # This is equivalent to the voronoi tesselation described in Krakow, + # but done in rodrigues space to take advantage of rectilinear planes. + dp = get_distinguished_points(start, end) + # These large cell normals are always one of the 15 from figure 5 of + # krakow2017 + large_cell_normals = _get_large_cell_normals(dp) + + # Step 3: (only for misorientations) restrict the domain to the + # fundamental sector of the pole figure of the shared symmetries. + disjoint = start & end fz = disjoint.fundamental_zone() fz_normals = Rotation.from_axes_angles(fz, np.pi) normals = Rotation(np.concatenate([large_cell_normals.data, fz_normals.data])) diff --git a/orix/quaternion/quaternion.py b/orix/quaternion/quaternion.py index 84bbf638..f0e50485 100644 --- a/orix/quaternion/quaternion.py +++ b/orix/quaternion/quaternion.py @@ -1114,8 +1114,15 @@ def dot_outer(self, other: Quaternion) -> np.ndarray: dots = np.tensordot(self.data, other.data, axes=(-1, -1)) return dots - def mean(self) -> Quaternion: - """Return the mean quaternion with unitary weights. + def mean(self, weights: np.ndarray | None = None) -> Quaternion: + r"""Return the mean quaternion. + + Parameters + ---------- + weights + An optional array of weights for calculating a weighted + average instead of the unweighted mean. Must be the same + size as the quaternion array. Returns ------- @@ -1124,14 +1131,41 @@ def mean(self) -> Quaternion: Notes ----- - The method used here corresponds to Equation (13) in - https://arc.aiaa.org/doi/pdf/10.2514/1.28949. + The method used here corresponds to equations 12 and 13 in + :cite:`markley_averaging_2007`, which simplifies to the + following two equations: + + .. math:: + + M = \sum_{n=1}^{n}{wq_i \cdot q_i^T} + + q_{mean} = \max_{q \in SO(3)} \left( q^T \cdot M \cdot q \right) + + Which gives the Frobenius norm of rotation space (SO(3)). This + is different from the following Euclidean-style equation + occasionally used elsewhere in crystallography: + + .. math:: + + q_{mean} = norm\left( \frac{\sum_{n=1}^{n}{wq_i }}{n} \right) + + Which gives the Frobenius norm of quaternion space (SU(2)). + These are roughly identical for closely clustered quaternions, + but as explained in :cite:`markley_averaging_2007`, + increasingly deviate as the spread in the data increases. """ Q = self.flatten().data.T - QQ = Q.dot(Q.T) + if weights is not None: + weights = np.asanyarray(weights).flatten()[:, np.newaxis] + QQ = Q.dot(weights * Q.T) + else: + QQ = Q.dot(Q.T) w, v = np.linalg.eig(QQ) - w_max = np.argmax(w) - return self.__class__(v[:, w_max]) + v_mean = v[:, np.argmax(w)] + # flip if necessary. + if v_mean[0] < 0: + v_mean = v_mean * -1 + return self.__class__(v_mean) def outer( self, diff --git a/orix/quaternion/rotation.py b/orix/quaternion/rotation.py index 79dd9a0e..e5a4a719 100644 --- a/orix/quaternion/rotation.py +++ b/orix/quaternion/rotation.py @@ -143,6 +143,12 @@ def __getitem__(self, key) -> Rotation: R.improper = self.improper[key] return R + def __setitem__(self, key, value: np.ndarray) -> None: + if isinstance(value, Rotation): + self._data[key] = value._data + else: + self.data[key] = value.data + def __invert__(self) -> Rotation: R = super().__invert__() R.improper = self.improper diff --git a/orix/tests/test_crystal_map/test_crystal_map.py b/orix/tests/test_crystal_map/test_crystal_map.py index 22237334..4d5507a1 100644 --- a/orix/tests/test_crystal_map/test_crystal_map.py +++ b/orix/tests/test_crystal_map/test_crystal_map.py @@ -678,10 +678,10 @@ def test_orientations(self, crystal_map_input, phase_list): @pytest.mark.parametrize( "point_group, rotation, expected_orientation", [ - (C2, [(0.6088, 0, 0, 0.7934)], [(-0.7934, 0, 0, 0.6088)]), - (C3, [(0.6088, 0, 0, 0.7934)], [(-0.9914, 0, 0, 0.1305)]), - (C4, [(0.6088, 0, 0, 0.7934)], [(-0.9914, 0, 0, -0.1305)]), - (O, [(0.6088, 0, 0, 0.7934)], [(-0.9914, 0, 0, -0.1305)]), + (C2, [(0.6088, 0, 0, 0.7934)], [(0.7934, 0, 0, -0.6088)]), + (C3, [(0.6088, 0, 0, 0.7934)], [(0.9914, 0, 0, -0.1305)]), + (C4, [(0.6088, 0, 0, 0.7934)], [(0.9914, 0, 0, 0.1305)]), + (O, [(0.6088, 0, 0, 0.7934)], [(0.9914, 0, 0, 0.1305)]), ], ) def test_orientations_symmetry(self, point_group, rotation, expected_orientation): diff --git a/orix/tests/test_plot/test_stereographic_plot.py b/orix/tests/test_plot/test_stereographic_plot.py index e4ec91b7..2b6762d6 100644 --- a/orix/tests/test_plot/test_stereographic_plot.py +++ b/orix/tests/test_plot/test_stereographic_plot.py @@ -345,7 +345,7 @@ def test_format_coord(self): plt.close("all") - def test_empty_scatter(self): + def test_out_of_range_scatter(self): v = Vector3d([0, 0, 1]) _, ax = plt.subplots(subplot_kw=dict(projection=PROJ_NAME)) @@ -358,6 +358,22 @@ def test_empty_scatter(self): plt.close("all") + def test_empty_scatter(self): + v = Vector3d.empty() + + _, ax = plt.subplots(subplot_kw=dict(projection=PROJ_NAME)) + ax.hemisphere = "lower" + + # Not plotted since the vector is empty + ax.scatter(v) + ax.text(v, s="1") + assert len(ax.texts) == 0 + + # Check a stereographic plot was still made + assert ax._has_collection("_stereographic_polar_grid")[0] + + plt.close("all") + @pytest.mark.parametrize("shape", [(5, 10), (2, 3)]) def test_multidimensional_vector(self, shape): n = np.prod(shape) diff --git a/orix/tests/test_quaternion/test_misorientation.py b/orix/tests/test_quaternion/test_misorientation.py index ab2a1060..d1761ec5 100644 --- a/orix/tests/test_quaternion/test_misorientation.py +++ b/orix/tests/test_quaternion/test_misorientation.py @@ -21,13 +21,14 @@ import numpy as np import pytest from scipy.spatial.transform import Rotation as SciPyRotation +from scipy.stats import norm from orix._utils.constants import VisibleDeprecationWarning from orix.crystal_map import Phase -from orix.quaternion import Misorientation +from orix.quaternion import Misorientation, Quaternion # isort: off -from orix.quaternion.symmetry import C1, D6, Oh, _groups +from orix.quaternion.symmetry import C1, C2h, C3, C2v, D6, T, Oh, _groups # isort: on from orix.vector import Miller, Vector3d @@ -186,15 +187,36 @@ def test_random(self): M3 = Misorientation.random(symmetry=(Oh, D6)) assert M3.symmetry == (Oh, D6) - # TODO: Remove after v0.15.0 is released - def test_map_into_symmetry_reduced_zone_deprecation_message(self): - M = Misorientation.random() - M.symmetry = (Oh, Oh) - with pytest.warns( - VisibleDeprecationWarning, - match=( - r"Function `map_into_symmetry_reduced_zone\(\)` is deprecated and will " - r"be removed in version 0.15. Use `reduce\(\)` instead." - ), - ): - _ = M.map_into_symmetry_reduced_zone() + def test_mean(self): + # create a random loosely clustered group of misorientations + np.random.seed(2319) + qu_data = np.stack([norm.rvs(i, 0.2, 20) for i in [0.3, 0.1, 0.2, 0.3]]).T + # The symmetries tested are Identity, a Laue, an improper, and an + # improper with no inversion point, which tests all combinations of + # OrientationRegion.from_symmetry() if/then logic. + syms = [C1, C2h, C3, C2v] + for start in syms: + for end in syms: + m = Misorientation(qu_data, symmetry=(start, end)).reduce() + # Test every variant of inputs works + rough = m.reduce().mean(ignore_symmetry=True) + p_mean, p_neigh = m.mean(include_improper=False, return_neighbors=True) + f_mean, f_neigh = m.mean( + include_improper=True, return_neighbors=True, verbose=True + ) + # for the three above calls, the deviation of the mean might + # lessen as symmetry constrains are added, and will never + # increase + r_dp = np.mean(np.abs(Quaternion(rough.data).dot(Quaternion(m.data)))) + f_dp = np.mean( + np.abs(Quaternion(f_mean.data).dot(Quaternion(f_neigh.data))) + ) + p_dp = np.mean( + np.abs(Quaternion(p_mean.data).dot(Quaternion(p_neigh.data))) + ) + assert r_dp <= p_dp + assert p_dp <= f_dp + # Test weighting + m1 = m[[0, 0, 0, 1, 2, 2, 4]] + m2 = m[:5] + m1.mean() == m2.mean(weights=[3, 1, 2, 0, 1]) diff --git a/orix/tests/test_quaternion/test_orientation.py b/orix/tests/test_quaternion/test_orientation.py index 110b5199..3d2e85f8 100644 --- a/orix/tests/test_quaternion/test_orientation.py +++ b/orix/tests/test_quaternion/test_orientation.py @@ -89,13 +89,13 @@ def test_quaternion_subclasses_copy_constructor_casting(): ([(1, 0, 0, 0)], T, [(1, 0, 0, 0)]), ([(1, 0, 0, 0)], O, [(1, 0, 0, 0)]), # 7pi/12 -C2-> # 7pi/12 - ([(0.6088, 0, 0, 0.7934)], C2, [(-0.7934, 0, 0, 0.6088)]), + ([(0.6088, 0, 0, 0.7934)], C2, [(0.7934, 0, 0, -0.6088)]), # 7pi/12 -C3-> # 7pi/12 - ([(0.6088, 0, 0, 0.7934)], C3, [(-0.9914, 0, 0, 0.1305)]), + ([(0.6088, 0, 0, 0.7934)], C3, [(0.9914, 0, 0, -0.1305)]), # 7pi/12 -C4-> # pi/12 - ([(0.6088, 0, 0, 0.7934)], C4, [(-0.9914, 0, 0, -0.1305)]), + ([(0.6088, 0, 0, 0.7934)], C4, [(0.9914, 0, 0, 0.1305)]), # 7pi/12 -O-> # pi/12 - ([(0.6088, 0, 0, 0.7934)], O, [(-0.9914, 0, 0, -0.1305)]), + ([(0.6088, 0, 0, 0.7934)], O, [(0.9914, 0, 0, 0.1305)]), ], indirect=["orientation"], ) @@ -393,9 +393,7 @@ def test_from_matrix_symmetry(self): assert o1.symmetry.name == "1" o2 = Orientation.from_matrix(om, symmetry=Oh) o2 = o2.reduce() - assert np.allclose( - o2.data, np.array([1, 0, 0, 0] * 2 + [-1, 0, 0, 0] * 2).reshape(4, 4) - ) + assert np.allclose(o2.data, np.array([1, 0, 0, 0] * 4).reshape(4, 4)) assert o2.symmetry.name == "m-3m" o3 = Orientation(o1.data, symmetry=Oh) o3 = o3.reduce() @@ -727,8 +725,20 @@ def test_reduce_verbose(self): o2 = o.reduce(verbose=True) assert np.allclose(o1.data, o2.data) - @pytest.mark.flaky(reruns=3) + # @pytest.mark.flaky(reruns=3) def test_reduce_all_groups(self): + np.random.seed(2319) for group in _groups: ori = Orientation.random(symmetry=group) + # NOTE: in the past, this test has needed a flaky flag. PR 669 + # should have fixed the underlying issue, but leaving a note here + # in case a similar bug occurs in the future/ assert np.isclose(ori.angle_with(ori.reduce()), 0) + + def test_dot(self): + # assert shape is correct + o1 = Orientation.random(shape=(2, 3, 4), symmetry=T) + o2 = Orientation.random(shape=(5, 6), symmetry=D3) + o12 = o2.dot_outer(o1) + assert o12.shape == (2, 3, 4, 5, 6) + assert np.all(o12[1, 2, :, 3, 2] == o2[3, 2].dot(o1[1, 2, :])) diff --git a/orix/tests/test_quaternion/test_orientation_region.py b/orix/tests/test_quaternion/test_orientation_region.py index bbe63b56..d7fa2c00 100644 --- a/orix/tests/test_quaternion/test_orientation_region.py +++ b/orix/tests/test_quaternion/test_orientation_region.py @@ -121,6 +121,9 @@ def test_get_distinguished_points(s1, s2, expected): def test_get_large_cell_normals(s1, s2, expected): n = _get_large_cell_normals(s1, s2) assert np.allclose(n.data, expected, atol=1e-3) + dp = s1.outer(s2).antipodal.unique(antipodal=False) + n_from_dp = _get_large_cell_normals(dp=dp) + assert np.allclose(n_from_dp.data, expected, atol=1e-3) def test_coverage_on_faces(): @@ -139,13 +142,8 @@ def test_coverage_on_faces(): (Csz, Ci), (C1, C1), (Ci, Ci), + (Csz, Csz), ], ) def test_get_proper_point_groups(Gl, Gr): _ = get_proper_groups(Gl, Gr) - - -def test_get_proper_point_group_not_implemented(): - """Double inversion case not yet implemented.""" - with pytest.raises(NotImplementedError): - _ = get_proper_groups(Csz, Csz) diff --git a/orix/tests/test_quaternion/test_quaternion.py b/orix/tests/test_quaternion/test_quaternion.py index 92eab8f9..d9e8e8b2 100644 --- a/orix/tests/test_quaternion/test_quaternion.py +++ b/orix/tests/test_quaternion/test_quaternion.py @@ -468,7 +468,8 @@ def test_from_to_matrix(self): def test_from_euler_to_matrix_from_matrix(self, eu): q = Quaternion.from_euler(eu.reshape(5, 2, 3)) - assert np.allclose(Quaternion.from_matrix(q.to_matrix()).data, q.data) + q_from_to = Quaternion.from_matrix(q.to_matrix()) + assert np.allclose(q_from_to.data, q.data, atol=1e-6) def test_from_matrix_to_euler_from_euler_to_matrix(self, eu): rot_180x = np.diag([1, -1, -1]) diff --git a/orix/tests/test_vector3d/test_vector3d.py b/orix/tests/test_vector3d/test_vector3d.py index 2b20dfa0..d6250606 100644 --- a/orix/tests/test_vector3d/test_vector3d.py +++ b/orix/tests/test_vector3d/test_vector3d.py @@ -496,6 +496,15 @@ def test_get_nearest(): v.get_nearest(v_ref) +def test_tuples(): + v = Vector3d(np.zeros([10, 3])) + tup = v._tuples + assert len(tup) == 1 + v = Vector3d(np.arange(30).reshape([10, 3])) + tup = v._tuples + assert len(tup) == 10 + + class TestSpareNotImplemented: def test_radd_notimplemented(self, vector): with pytest.raises(TypeError):