Behavior at singularities

[Fab7c4]

I do not want to go into discussions of Euler angles especially at singularities
Anyways we did look into the behavior at the singularity and compared it to different implementations

And for a minimal example:

import rowan
import numpy as np

def test_at_singularity():
    rpy  = np.array((np.radians(90.0), np.radians(90.0), np.radians(0.0)))
    ypr  = np.array([rpy[-1], rpy[1], rpy[0]])

    keywords_list = [{"convention": "xyz", "axis_type": "extrinsic"},
                     {"convention": "zyx", "axis_type": "intrinsic"},]
    
    values = [rpy, ypr]

def test_at_singularity():
    rpy  = np.array((np.radians(90.0), np.radians(90.0), np.radians(0.0)))
    ypr  = np.array([rpy[-1], rpy[1], rpy[0]])

    keywords_list = [{"convention": "xyz", "axis_type": "extrinsic"},
                     {"convention": "zyx", "axis_type": "intrinsic"},]
    
    values = [rpy, ypr]

    for euler, keywords in zip(values, keywords_list):
        print(f'Convention: {keywords["convention"]}, Axis type: {keywords["axis_type"]}')
        print("Euler", np.degrees(euler))
        quat = rowan.from_euler(*euler, **keywords)
        print("Quat", quat)
        euler_back = rowan.to_euler(quat, **keywords)
        print("Euler back", np.degrees(euler_back))
        quat_back = rowan.from_euler(*euler_back, **keywords)
        print("Quat back", quat_back)
        print("Symmetric distance", rowan.geometry.sym_intrinsic_distance(quat, quat_back))
        print("")
              

if __name__ == "__main__":
    test_at_singularity()

We obtain:

Convention: xyz, Axis type: extrinsic
Euler [90. 90.  0.]
Quat [ 0.5  0.5  0.5 -0.5]
Debug: at singularity
Euler back [-90.  90.   0.]
Quat back [ 0.5 -0.5  0.5  0.5]
Symmetric distance 1.5707963267948963

Convention: zyx, Axis type: intrinsic
Euler [ 0. 90. 90.]
Quat [ 0.5  0.5  0.5 -0.5]
Debug: at singularity
Euler back [-90.  90.   0.]
Quat back [ 0.5  0.5  0.5 -0.5]
Symmetric distance 0.0

An expected result for the intrinsic case yet not for the extrinsic case.
Now we saw if you adapt the in the case also the use of the multipliers as you do for 'xyz' intrinsic convention.
Things seem to work as expected:

    elif convention == "zyx":
        beta = np.arcsin(-mats[..., 2, 0])
        multiplier = mats[..., 2, 0] if axis_type == "extrinsic" else 1
        where_zero = np.isclose(np.cos(beta), 0, atol=atol)
        gamma = np.where(where_zero, 0, np.arctan2(mats[..., 2, 1], mats[..., 2, 2]))
        alpha = np.where(where_zero, 0, np.arctan2(mats[..., 1, 0], mats[..., 0, 0]))
        zero_terms = np.arctan2(multiplier *-mats[..., 0, 1], mats[..., 1, 1])

Convention: xyz, Axis type: extrinsic
Euler [90. 90.  0.]
Quat [ 0.5  0.5  0.5 -0.5]
Debug: at singularity
Euler back [90. 90.  0.]
Quat back [ 0.5  0.5  0.5 -0.5]
Symmetric distance 0.0

Convention: zyx, Axis type: intrinsic
Euler [ 0. 90. 90.]
Quat [ 0.5  0.5  0.5 -0.5]
Debug: at singularity
Euler back [-90.  90.   0.]
Quat back [ 0.5  0.5  0.5 -0.5]
Symmetric distance 0.0

https://github.com/glotzerlab/rowan/blob/7561ec80f4f970f8dd60b33fe9cbb002044b1a7c/rowan/functions.py#L859C32-L859C32