Orientation from align vectors

[anderscmathisen]

I have been trying to find a way to initialize an Orix Orientation object from an indexed diffraction pattern. My thinking is that if we know that u1 || v1 and u2 || v2 where u are Crystal directions (from the indexed diffraction pattern) and v are the sample directions (known from the x and y direction of the image of the diffraction pattern), then enough information is provided to determine the orientation of the crystal.

Scipys rotation object has a method for initializing a Rotation object from two sets of two vectors, where the rotation is described as the rotation which takes one from the first set of two vectors, to the second. It does this using Kabsch algorithm (https://en.wikipedia.org/wiki/Kabsch_algorithm).

Essentially what I am looking for is a similar function in Orix. I have made a very simple wrapper for scipys function to get a orix Orientation object:

def Orix_rotation_from_align_vectors(
    sample_vectors: Vector3d, crystal_millers: Vector3d, symmetry = None
) -> Orientation:
    vec1 = sample_vectors.unit.data
    vec2 = crystal_millers.unit.data
    scipy_rotation, error = ScipyRotation.align_vectors(vec2, vec1)
    rotation_matrix = scipy_rotation.as_matrix()
    
    return Orientation.from_matrix(rotation_matrix, symmetry)

I have also tried to code up some algorithm to get this Orientation myself. If the vectors in both sets we are trying to align are normal to each other, then an algorithm which first finds the rotation g1 which takes u1 to v1, and then finds the rotation g2 which takes the vector (g1 * u1) to v2 (which will just be a rotation about the vector v1) should work as the full rotation will be g = g2 * g1. However, if the vectors in the two sets are not normal, it seems to me that this algorithm will not work.

Would the wrapper above be a useful addition to the Orientation object in Orix? Or do you have any suggestions to how such a function could be implemented in Orix in the best way possible?

[hakonanes]

Hi @anderscmathisen,

thank you for suggesting this functionality. Pretty sure there are others in your situation, so such a class method should be useful to have in orix! If others agree, I suggest you open a PR adding this method to the Orientation class as a class method. I can help you with tests, documentation, formatting etc. (although everything you need should be listed in the contributing guide).

if the vectors in the two sets are not normal, it seems to me that this algorithm will not work.

I believe you are right there. I suggest to just wrap scipy.spatial.transform.Rotation.align_vectors() (docs) in this case. I know both @harripj and @argerlt have used SciPy's rotation functionality before and might be familiar with this method (I'm not).

Regarding the implementation, we could return the other returns from SciPy if the user asked for it in a return_* boolean parameter, and also allow passing the other parameters in align_vectors() as keyword arguments? And how about calling the parameters a and b, since we don't really put a restriction on what the sample_vectors and crystal_vectors are, right? Of course, examples in the documentation of using the method can show the use case you're mentioning. What do you think?

[anderscmathisen]

Hi,

Yes, I definetively agree that we should allow users to send additional possible args / kwargs to the scipy.spatial.transform.Rotation.align_vectors() function, as well as returning values such as the error upon request.

It also makes sense to be more general with the parameter names than the specific case of crystallography.

[argerlt]

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.

Two suggestions, take them or leave them:

1. I would suggest adding this functionality to the Quaternion class as a Quaternion.from_vectors or Quaternion.from_aligned_vectors function, so it can be used to define a rotation, quaternion, orientation, and misorientation via inheritance. Would also align better with Drop support for Python 3.6 #343, Purpose for neo_eulerian classes #374, and Selection of data from cluster analysis #375

2. Just an FYI in case you haven't fallen into this trap before, There is no universally accepted way to map rotations to quaternions. You can pick one of 8 equally valid hemispheres of the 4d unit quaternion hypersphere, 6 permutations of xyz on each hemisphere, and choose active or passive representation, giving 96 equally valid ways to map rotations to quaternion space. Orix and scipy use different mappings, and more confusingly, have different mappings from quaternions to rotation matrices as well. Case and point, try running the following code and see how they pop up differently:

from orix.quaternion import Orientation
from scipy.spatial.transform inport Rotation as R

euler = np.array([15,32,41])*np.pi/180

orix_quat = Orientation.from_euler(e)
scipy_quat = R.from_euler('zyz',e).as_quat()
print(orix_quat)
print(scipy_quat)

notably, orix follows the (s, v1,v2,v3) convention while scipy orders theirs as (v1,v2,v3,s). These are two of those 8 valid hemispheres. I THINK this is the magic formula to convert between the two, but no guarantees:
correct_scipy_quat = R.from_euler('zxz',e).as_quat()[(3,0,1,2),]*np.array([1,-1,1,-1])

With that in mind, I would instead suggest the following version of your function:

def Orix_rotation_from_align_vectors(vec1,vec2):
    vec1 = sample_vectors.unit.data
    vec2 = crystal_millers.unit.data
    scipy_rotation, error = ScipyRotation.align_vectors(vec2, vec1)
    scipy_quat = scipy_rotation.as_quat    
    # insert math here to change from scipy to orix quaternion. 

    return Quaternion(orix_quat )
'''

Since sicpy rotations are stored as quaternions internally, this lets you avoid the added confusion from transferring from scipy quaternion to scipy rotation matrix to orix rotation matrix to orix quaternion.  

[anderscmathisen]

Hi,

Yes, I agree that this function is perhaps best implemented in the Quaternion class and inherited to the others. I will give the conversion calculations a try :)

[anderscmathisen]

I am currently working on a PR for this.

First of all, I have started to implement the function as part of the Rotation class and not Quaternion. This might be me not understanding the philosophy of Orix, but I noticed that the Quaternion class defined in orix/quaterion/quaternion.py did not have any similar functions such as from_euler() ect, and so it felt more natural to add the function along with the other from_...() functions of Rotation class.

Second, most of the from_...() functions such as from_euler() support creating arrays of multiple orientations. That is, if you do Rotation.from_euler([[0,0,np.pi], [0,np.pi,0]]) you get an object with two rotations. The function scipy.spatial.transform.Rotation.align_vectors() does not really support this type of behavior. The vector sets a and b can contain an arbitrary amount of vectors, but only a single rotation going from vectors in b to a will be retuned. With some more wrapping, it should however be possible to implement the orix function to return multiple rotations.

One possible implementation is to have a and b have a shape of (M,N,3), where M is the number of rotations we want returned and N is the number of vectors in each set.

rotations = []
 for a_m, b_m in a, b:
    rotation = scipy.align_vectors(a_m,b_m)
    orix_rotation = convert_to_orix(rotation)
    rotations.append(orix_rotation)

return rotations

Do you have any suggestions here? Is it important to support creating arrays of rotations?

[hakonanes]

Is it important to support creating arrays of rotations?

I'd say no, since your use case is for a single orientation. I think best practice in terms of implementation and maintainability is to implement what you need, and not what you anticipate you or someone else will need in the future.

This might be me not understanding the philosophy of Orix

In general, if there is no particular reason to not add it to the parent class, it should go in the parent class. Rotation adds the notion of an improper rotation which Quaternion does not (see the docs).

We should wrap Quaternion.from_align_vectors() in Orientation.from_align_vectors() accepting symmetry, and Misorientation.from_align_vectors() accepting a symmetry tuple, but this can be discussed in the coming PR.

I actually think Quaternion.from_scipy_rotation() could be useful as its own class method: you could then call this within the from_align_vectors() method. @argerlt might have input to both of my last two points here.