Add Quaternion.to_scipy_rotation?

[ericpre]

I am wondering if it would be sensible to add the functionality to convert orix quaternion to scipy rotation (the reverse is already possible with Rotation.from_scipy_rotation). I would find it useful to define rotation to be used in the rotation of block waves in abtem.

from orix.crystal_map import Phase
from orix.quaternion import Rotation
from orix.vector import Miller

phase = Phase(point_group="mmm")
t = Miller.from_highest_indices(phase, uvw=[1, 1, 1])
t = t.in_fundamental_sector()
t = t.unit.unique(use_symmetry=True).round()
print(t)

vz = Miller(uvw=[0, 0, 1], phase=t.phase)
R = Rotation.identity(t.size)
for i, t_i in enumerate(t):
    R[i] = Rotation.from_align_vectors(t_i, vz)
print(R)

from scipy.spatial.transform import Rotation as ScipyRotation
from orix.quaternion import Rotation as Rotation

scipy_R = ScipyRotation.from_euler('ZXZ', R.to_euler()).as_euler("xyz")
print("Euler angle (Scipy rotation):")
print(scipy_R)

[hakonanes]

This is a very good idea. The use of the intrinsic letters "ZXZ" (rotation of reference frame, not object) seems correct to me.

Do you have time to implement it?

As far as I can tell, sub-classes of Quaternion don't have to overwrite this method as they will just use Quaternion.to_euler().

A test must ensure a consistent conversion loop.

And a nice example in the gallery would be great!

[ericpre]

I am not familiar with the subtilties of the rotation business and I think that will take me some time to figure out what are the various convention and how they are used by scipy and orix. I suspect that it would be easier for someone that understand better all this, in order word, it is fairly unlikely that I will do it.

[hakonanes]

Sure, I can do it. Most likely within a 1,5 weeks.

[argerlt]

I think there is a faster and less convoluted way to do this where you just pass the reordered quaternion. For testing, the easiest approach would probably be passing the quaternion unittest data through Rotation.to_scipy_rotation().as_rotvec() and comparing to `qu2ax'.

I can't do proper pull requests till the new year, but if this is still open in January, I'll add it.

As a side note, converting via Euler angle triplets can be a headache. Gimble lock and variable accuracy aside, switching passive to active also switches the ordering of the triplets. There are also a lot of historical terminology problems (for example, what orix and MTEX call bunge angles aren't the triplet originally used by H.J. bunge, see section 1.8.2 of C-S Man's book for details. Using the code below though, the equivalent angle triplets are:

orix_r = Rotation.from_axes_angles([1,2,3],np.pi/9)
scipy_r = SciPyRotation(orix_r.data[:,(1,2,3,0)]*[1,1,1,-1])

orix_r.to_euler()
scipy_r.as_euler('ZXZ')

Out[185]: array([[-2.17488927,  0.20792361,  1.8939986 ]])

Out[186]: array([[4.10829604, 0.20792361, 1.8939986 ]])

note even here there is confusion, as scipy defaults to the smallest absolute angle magnitude whereas orix chooses positive angles only.

def to_scipy_rotation(self,invert_to_active: bool=False):
    """ Convert a passive scalar-vector orix transformation into an active
    vector-scalar scipy transfomation

    IMPORTANT CONVERSION CONVENTION INFORMATION    
    Assume the following notation:
        w = cos(theta / 2)
        x = sin(theta / 2) * n_x
        y = sin(theta / 2) * n_y
        z = sin(theta / 2) * n_z
    There are at least 8 equally valid ways to map a rigid body rotation to
    the quaternion sphere. of those, orix uses a passive scalar-vector order
    [w, x, y, z], describing the rotation of a reference frame into an object.
    Scipy, by default, uses an active vector-scalar [x, y, z, w], describing
    the rotation of an object to a reference frame.
    
    by default, `to_scipy_rotation' will pass [x, y, z, w] into scipy, which
    technically inverts the rotation, but gives the expected results when
    comparing rotation vectors and euler angles:

        import scipy.spatial.transform.rotation as SciPyRotation
        import orix.quaternion.Rotation as Rotation
        orix_r = Rotation.from_axes_angles([1,2,3],np.pi/9)
        scipy_r = orix_r.to_scipy_rotation()
        # axes
        scipy_r.as_rotvec()/(scipy_r.as_rotvec().data[0,0])
        orix_r.axis/orix_r.axis.x

    >>> array([[1., 2., 3.]])
    >>> Vector3d (1,)
        [[1. 2. 3.]

         # angles
        orix_r.angle*180/np.pi
        scipy_r.magnitude()*180/np.pi

    >>> array([20.]
    >>> array([20.]

      However, since scipy is ACTIVE and orix is PASSIVE, the 'equivalent'
    representation is to pass `invert_to_active=True'. this gives identical
    axes but inverted angles, which can be confusing to users unfamiliar with
    active/passive conventions.
        '''    
        scipy_r = orix_r.to_scipy_rotation(True)
        # axes
        scipy_r.as_rotvec()/(scipy_r.as_rotvec().data[0,0])

    >>> array([[1., 2., 3.]])

         # angles
        scipy_r.magnitude()*180/np.pi

    >>> array([-20.]
    '''
    
    Parameters
    ----------
    invert_to_active : bool, optional
        If True, inverts passive (intrinsic) orix quaternions to active 
        (extrinsic) scipy quaternions. If False, no conversion is made and
        data is passed as-is, which technically inverts the rotation angle.
        The default is False.

    Returns
    -------
    scipy.spatial.transform.Rotation
        scipy rotation object.

    """
    # scipy uses the vector-scalar quaternion instead of the scalar-vector.
    scipy_quat = self.data[:,(1,2,3,0)]
    # scipy is active, orix is passive. alllow the option to invert as
    # desired.
    if invert_to_active:
        scipy_quat = scipy_quat* [1,1,1,-1]
    return SciPyRotation.from_quat(scipy_quat)
    ```