Add class method from_align_vectors() to Quaternion, Misorientation and Orientation

[anderscmathisen]

Description of the change

This PR is a result of an ongoing recent discussion.

I have made functions to wrap scipy.spacial.transform.Rotation.align_vectors() to initialize Quaternion, Misorientation or Orientation objects which rotates a set of vectors b into another set a. For the Misorientation and Orientation objects, symmetry arguments are also accepted.

For the Quaternion class, I have also made a function which initializes a Quaternion object from a Scipy Rotation object.

The implementation of multiple return values in my code is not very elegant, and suggestions for making it better is appreciated.

In the discussion, the conversion from Scipy quaternions to Orix quaternions is discussed. I did look a little into this, and got very confused, but the current conversion in Quaternion.from_scipy_rotation() seems to work in the sense that the resulting orientation from the function Quaternion.from_align_vectors() indeed transforms vectors from b to a.

Progress of the PR

• [X (ish) ] Docstrings for all functions
• Unit tests with pytest for all lines
• Clean code style by running black via pre-commit

Minimal example of the bug fix or new feature

>>> from orix.quaternion import Orientation
>>> from orix.vector import Vector3d, Miller
>>>
>>> crystal_millers = Miller([[2,-1,0],[0,0,1]])
>>> sample_vectors = Vector3d([[3,1,0],[-1,3,0]])
>>>
>>> ori = Orientation.from_align_vectors(crystal_millers, sample_vectors)
>>> print(crystal_millers)
>>> print(ori * sample_vectors.unit * crystal_millers.length)

Output:

Miller (2,), point group None, xyz
[[ 2 -1  0]
 [ 0  0  1]]
Vector3d (2,)
[[ 2. -1. -0.]
 [-0. -0.  1.]]

For reviewers

• The PR title is short, concise, and will make sense 1 year later.
• [n/a] New functions are imported in corresponding __init__.py.
• New features, API changes, and deprecations are mentioned in the unreleased
  section in CHANGELOG.rst.
• Contributor(s) are listed correctly in __credits__ in orix/__init__.py and in
  .zenodo.json.

[hakonanes]

Thank you for taking the time to set up this PR, @anderscmathisen!

This first review is just of the quaternion methods. Will come back to the others once I've found some tests for them to convince myself of their correctness. Then we can write some tests and add some examples to the docstrings.

I find your implementation in good shape. I have many smaller suggestions for changes to keep the code style of the rest of the orix code base (see the code style guide) and some clarifications of the docstrings. Feel free to question my comments!

[hakonanes]

Thank you for resolving most of my comments. Below are a few more comments.

I'm confident that Orientation.from_align_vectors() is OK, based on the following test, which we could simplify and add to the test suite (as multiple smaller unit tests):

from diffpy.structure import Lattice, Structure
import numpy as np
from orix.crystal_map import Phase
from orix.quaternion import Orientation, symmetry
from orix.vector import Miller, Vector3d


phases = {
    "cubic": Phase(point_group="m-3m"),
    "hexagonal": Phase(
        point_group="6/mmm",
        structure=Structure(lattice=Lattice(1, 1, 2, 90, 90, 120))
    ),
    "trigonal": Phase(
        point_group="-3m",
        structure=Structure(lattice=Lattice(4.9, 4.9, 5.4, 90, 90, 120))
    ),
    "tetragonal": Phase(
        point_group="4/mmm",
        structure=Structure(lattice=Lattice(0.5, 0.5, 1, 90, 90, 90))
    ),
    "orthorhombic": Phase(
        point_group="mmm",
        structure=Structure(lattice=Lattice(1, 2, 3, 90, 90, 90))
    ),
    "monoclinic": Phase(
        point_group="2/m",
        structure=Structure(lattice=Lattice(1, 2, 3, 90, 93, 90))
    ),
    "triclinic": Phase(
        point_group="-1",
        structure=Structure(lattice=Lattice(1, 2, 3, 70, 100, 120))
    )
}
phase = phases["triclinic"]

# Define a reference orientation (goal)
ori_ref = Orientation.from_axes_angles([1, 1, 1], np.pi / 4, phase.point_group)

# Specify two crystal directions (any will do)
uvw = [[2, 1, 1], [1, 3, 1]]
vec_c = Miller(uvw=uvw, phase=phase)

# Find out where these directions in the reference orientation (crystal)
# point in the sample reference frame
vec_r = Vector3d(~ori_ref * vec_c)

# Plot the reference orientation directions as empty circles
fig = vec_r.scatter(
    ec=["r", "b"],
    s=100,
    fc="none",
    grid=True,
    axes_labels=["X", "Y"],
    return_figure=True,
    hemisphere="both",
)

# Add some randomness to the sample directions (0 error magnitude gives
# exact result)
mag_err = 0
vec_err = Vector3d(np.random.normal(0, mag_err, 3))
vec_r_err = vec_r + vec_err
angle_err = np.rad2deg(vec_r_err.angle_with(vec_r)).mean()
print("Vector deviation: ", angle_err)

# Obtain the rotation which aligns the crystal directions with the
# sample directions (both sample to crystal, r2c, and vice versa)
ori_new_r2c, err = Orientation.from_align_vectors(
    vec_c, vec_r_err, return_rmsd=True
)
ori_newc2r = ~ori_new_r2c
print("Error: ", err)

# Assert equality
print(np.allclose(ori_ref.data, ori_new_r2c.data))
print(np.allclose(ori_new_r2c.data, (~ori_new_c2r).data))

# Plot the crystal directions in the new orientation
vec_r2 = Vector3d(~ori_new_r2c * vec_c)
vec_r2.scatter(c=["r", "b"], figure=fig)

Output:

Vector deviation:  0.0
Error:  0.0
True
True

[anderscmathisen]

Thank you for the excellent feedback, @hakonanes. Glad to see the test you wrote works:)

Regarding the sorting of the parameters, I have inserted the symmetry argument in method from_align_vectors() for Orientation and Misorientation before the weights argument. The reason for this is that I, and I suspect most people, will use the symmetry argument much more than the weights argument.

[hakonanes]

Thanks again for accomodating my comments.

After having looked into (mis)orientation methods, I think we should narrow down their purpose so that they do what we intend them to do and nothing else. I've made two comments to that effect.

Here's another test I've made to convince myself of the use and correctness of Misorientation.from_align_vectors():

from diffpy.structure import Lattice, Structure
import numpy as np
from orix.crystal_map import Phase
from orix.quaternion import Misorientation, Orientation
from orix.vector import Miller, Vector3d


phases = {
    "cubic": Phase(point_group="m-3m"),
    "hexagonal": Phase(
        point_group="6/mmm",
        structure=Structure(lattice=Lattice(1, 1, 2, 90, 90, 120))
    ),
    "trigonal": Phase(
        point_group="-3m",
        structure=Structure(lattice=Lattice(4.9, 4.9, 5.4, 90, 90, 120))
    ),
    "tetragonal": Phase(
        point_group="4/mmm",
        structure=Structure(lattice=Lattice(0.5, 0.5, 1, 90, 90, 90))
    ),
    "orthorhombic": Phase(
        point_group="mmm",
        structure=Structure(lattice=Lattice(1, 2, 3, 90, 90, 90))
    ),
    "monoclinic": Phase(
        point_group="2/m",
        structure=Structure(lattice=Lattice(1, 2, 3, 90, 93, 90))
    ),
    "triclinic": Phase(
        point_group="-1",
        structure=Structure(lattice=Lattice(1, 2, 3, 70, 100, 120))
    )
}

# Specify two phases, possibly different
phase1 = phases["cubic"]
phase2 = phases["hexagonal"]

# Specify one orientation per crystal (phase)
ori_ref1 = Orientation.from_axes_angles(
    [1, 1, 1], np.pi / 3, symmetry=phase1.point_group
)
ori_ref2 = Orientation.from_axes_angles(
    [1, 3, 2], np.pi / 2, symmetry=phase2.point_group
)

# Get reference misorientation (goal). Transformations are composed from
# the right, so: crystal 1 -> sample -> crystal 2
mori_ref = Misorientation(
    ori_ref2 * (~ori_ref1), symmetry=(ori_ref1.symmetry, ori_ref2.symmetry)
)

# Specify two directions in the first crystal
vec_c1 = Miller(uvw=[[1, 1, 0], [0, 0, 1]], phase=phase1)

# Express the same directions with respect to the second crystal
vec_c2 = Miller(xyz=(mori_ref * vec_c1).data, phase=phase2)

# Add some randomness to the second crystal directions (0 error
# magnitude gives exact result)
mag_err = 0
vec_err = Miller(xyz=np.random.normal(0, mag_err, 3), phase=phase2)
vec_c2_err = Miller(xyz=(vec_c2 + vec_err).data, phase=phase2)
angle_err = np.rad2deg(vec_c2_err.angle_with(vec_c2)).mean()
print("Vector deviation: ", angle_err)

# Get the misorientation that aligns the directions in the first crystal
# with those in the second crystal
mori_new, err = Misorientation.from_align_vectors(vec_c2_err, vec_c1, return_rmsd=True)
print("Error: ", err)

# Plot the two directions in the (unrotated) first crystal as open
# circles
fig = vec_c1.scatter(
    ec=["r", "b"],
    s=100,
    fc="none",
    grid=True,
    axes_labels=["e1", "e2"],
    return_figure=True,
    hemisphere="both",
)

# Plot the two directions in the second crystal with respect to the
# first crystal's axes
(~mori_ref * vec_c2_err).scatter(c=["r", "b"], figure=fig)

# Check equality (False if an error is added)
print(np.allclose(mori_new.data, mori_ref.data))

Output:

Vector deviation:  0.0
Error:  0.0
True

[hakonanes]

@anderscmathisen, feel free to question my comments.

From #400 (comment):

This would be a super helpful function for me as well, as parent/child orientation relationships are often similarly given as vectors to avoid confusion, but MDF and ODF calculations require that same info as misorientations.

@argerlt, does the example I included above (#401 (review)) look like what you're after?

[argerlt]

@anderscmathisen, feel free to question my comments.

From #400 (comment):

This would be a super helpful function for me as well, as parent/child orientation relationships are often similarly given as vectors to avoid confusion, but MDF and ODF calculations require that same info as misorientations.

@argerlt, does the example I included above (#401 (review)) look like what you're after?

I believe so. I was actually coming here to give my suggestion for a from_scipy function, but it looks like the method @anderscmathisen worte passes your tests, so I'm happy!

@anderscmathisen, feel free to question my comments.

From #400 (comment):

This would be a super helpful function for me as well, as parent/child orientation relationships are often similarly given as vectors to avoid confusion, but MDF and ODF calculations require that same info as misorientations.

@argerlt, does the example I included above (#401 (review)) look like what you're after?

Yes it does. I was actually getting ready to submit a PR with my own from_scipy function, but @anderscmathisen 's version is is shorter and still passes all the suggested tests, so I have no feedback.

[hakonanes]

Yes it does.

Thank you for confirming.

[hakonanes]

I've got some comments for parameter naming and descriptions and some minor things for the implementation, otherwise I think it looks good. I suggest to add some docstring examples and tests after my comments are addressed.

[anderscmathisen]

I have added some examples in the docstrings of the methods. These are probably not the best examples, but probably work fine? I would however have liked to have a better method of displaying the results rather than np.allclose() and return value True, but the docstring examples test code gave me errors when I used print. Any ideas?

The test I added for Quaternion.from_scipy_rotation method fails for some reason, even though I think the code should be ok. This is probably related to what @argerlt wrote in the discussion, but Im not sure.

The other unit tests are probably a little too focused on checking the type of the return values, idk, these are my first tests. Feel free to give feedback. Perhaps it would have been better to just use the tests you have written in this PR, @hakonanes ?

[hakonanes]

Great that you've added docstring examples and tests! I find this PR close to completion.

[hakonanes]

Perhaps it would have been better to just use the tests you have written in this PR, @hakonanes ?

I think your tests are fine. I'd like to add my examples myself as examples in the doc examples section after this PR is ready, if that's OK with you.

[hakonanes]

Tests and examples are good! I have one final question which we must address, and after them I'm happy to make the final changes (add you to the credits list and update the changelog).

[hakonanes]

The failing test on Windows with Python 3.9 is handled in #402.

[hakonanes]

I'm happy with this PR now, great work @anderscmathisen!

Before merging, I'd like to wait for @harripj or @argerlt to give their take on my question at the bottom of #401 (comment).

[hakonanes]

How would you like your name listed in the credits @anderscmathisen, "Anders Christian Mathisen"? I'll add your name with NTNU as the affiliation in the Zenodo metadata file.

[hakonanes]

And this is the changelog entries I plan to add to your branch

- Creation of ``Quaternion``(s) (or instances of inheriting classes) from SciPy
  ``Rotation``(s).
- Creation of one ``Quaternion`` or ``Rotation`` by aligning sets of vectors in two
  reference frames, one ``Orientation`` by aligning sets of sample vectors and crystal
  vectors, and one ``Misorientation`` by aligning two sets of crystal vectors in two
  differente crystals.

[harripj]

Thanks for the ping @hakonanes, I will try to have a look a this before tomorrow!

[anderscmathisen]

How would you like your name listed in the credits @anderscmathisen, "Anders Christian Mathisen"? I'll add your name with NTNU as the affiliation in the Zenodo metadata file.

Anders Christian Mathisen is fine 👍🏼 The changelog also looks good to me.

[argerlt]

I'm happy with this PR now, great work @anderscmathisen!

Before merging, I'd like to wait for @harripj or @argerlt to give their take on my question at the bottom of #401 (comment).

I left a comment describing my approach to solving the testing problem. If you are interested and this issue is still open October 25th, I can write up a from_scipy function that uses rotations as the pass-between method, but I am swamped until then. Otherwise, @anderscmathisen or @hakonanes , feel free to write your own version based on my suggestion, or ignore my comment entirely in favor of a different approach.

If you DO want to write up scipy2orix via rot method yourself, I believe the method should be something like this:

@classmethod
def from_scipy_rotation(cls,assume_active_input=True):
    if assume_active_input == True:
        passive_matrix = np.linalg.inv(scipy_quat.as_matrix())
    else:
        passive_matrix =scipy_quat.as_matrix()
    orix_quat = Orientation.from_matrix(passive_matrix)

And the check would be something like:

# check bunge conversion
passive_bunges = (np.random.rand(1000,3) - 0.5)*np.array([np.pi*2,np.pi,np.pi*2])
active_bunges = passive_bunges[...,::-1]*(-1)
scipy_rot = R.from_euler('ZXZ',active_bunges).as_matrix()
orix_rot = Orientation.from_euler(passive_bunges).to_matrix()
assert(np.max(scipy_rot-orix_rot)<1E-10

but obviously much better formatted. Hope this helps!

[harripj]

Sorry for the delay all, and thanks for the interesting discussion! My only comment regarding inverting or not the Scipy rotation is that we document it in Notes of the docstring of from_scipy_rotation() so that anyone interested can see the interpretation (and then adjust for their use case as required).

As for deciding a default interpretation as discussed, given that SciPy rotations are active, I also think we should invert internally so that a SciPy rotation will be consistent in the orix framework.

[anderscmathisen]

I think my last commits here are consistent with our discussions on this PR. The current implementation uses matrix for scipy conversion and also inverts the scipy rotation. I added a short note in the docstring of from_scipy_rotation about the inversion. It's quite short, and perhaps not very easy to understand for new users?

[hakonanes]

After having run through my own examples/tests and sufficiently assured myself of the correctness of the use of from_align_vectors(), I'm happy with this PR. Good work, @anderscmathisen! And thanks to others for the discussion of conventions and the correct approach here.

I will make a follow-up PR adding the mentioned examples to the examples gallery. Feel free to do so yourself @anderscmathisen if you have a nice demonstration of the use of from_align_vectors() or any other demonstration!