Inconsistency with Orientation.angle_with

[soupmongoose]

I have encountered a potential bug with certain angles between two given orientations. The orientations start as euler angles and are then converted to orix Orientations.
Example:

oritest1 = Orientation.from_euler(np.radians([120.        ,  54.26282923,  45.        ]), symmetry=symmetry.Oh)
oritest2 = Orientation.from_euler(np.radians([ 0.        , 54.73561032, 45.        ]), symmetry=symmetry.Oh)

Then if I compare the results of

np.degrees((oritest1 - oritest2).angle.data)
np.degrees(oritest1.angle_with(oritest2).data)

I would expect the angle difference to be very small for both but only the first way gives the expected answer of 0.4727 while angle_with returns 30.0036. It seems to not take into account the 120degree symmetry in the first euler angle.

Further euler angle examples that give similarly strange results are:
[302. , -52.35178357 , -46.0036662 ] to [ 0. , 53.59948955 , 47.5 ]
[308.5 , -51.76150996 , -52. ] to [ 0. , 51.76150996 , 52. ]
[328.5 , -47.2657904 , -67.5 ] to [ 0. , 47.2657904 , 67.5 ]
These all have large differences when using angle_with and small, more expected, differences when using (Orientation - Orientation).angle.
However (Orientation - Orientation).angle also seems to have some flaws example:
[178.5 , 36.97748349 , 63.61592335] to [ 0. , 36.48968436 , 61.39659559]
Here Orientation.angle_with gives the expected result of 1.441 while (Orientation - Orientation).angle returns 39.268; seemingly not taking into account the 180degree symmetry in the first angle.

Hopefully this is enough to understand the issue, but if not I can provide more.

[hakonanes]

Thanks for reporting this @soupmongoose, it looks like a bug at a first glance.

Initialization of a Misorientation from two orientations:

orix/orix/quaternion/orientation.py

Lines 438 to 444 in 0fa14d1

	def __sub__(self, other):
	if isinstance(other, Orientation):
	# Call to Object3d.squeeze() doesn't carry over symmetry
	misorientation = Misorientation(self * ~other).squeeze()
	misorientation.symmetry = (self.symmetry, other.symmetry)
	return misorientation.map_into_symmetry_reduced_zone()
	return NotImplemented

Orientation.angle_with():

orix/orix/quaternion/orientation.py

Lines 539 to 553 in 0fa14d1

	def angle_with(self, other):
	"""The smallest symmetry reduced angle of rotation transforming
	this orientation to the other.
	Parameters
	----------
	other : orix.quaternion.Orientation
	Returns
	-------
	Scalar
	"""
	dot_products = self.unit.dot(other.unit).data
	angles = np.nan_to_num(np.arccos(2 * dot_products**2 - 1))
	return Scalar(angles)

I'll have a look into these routes and compare with results from MTEX when I find the time, hopefully within the weekend.

[harripj]

@soupmongoose thanks for the report and thanks @hakonanes ! I can shed some light on the angle_with situation, at the moment the function is written in a way that does not take into account the Orientation.symmetry, perhaps this should be updated!

For the final example I made some tests for your findings but was not able to reproduce the angles you quote in the final example:

>>> from orix.quaternion import Orientation
>>> import numpy as np
>>> o1 = Orientation.from_euler([178.5 , 36.97748349 , 63.61592335])
>>> o2 = Orientation.from_euler([ 0. , 36.48968436 , 61.39659559])
>>> print('C1 symmetry misori', np.rad2deg((o1 - o2).angle.data))
>>> print('C1 symmetry angle_with', np.rad2deg(o1.angle_with(o2).data))
C1 symmetry misori [136.89201375]
C1 symmetry angle_with [136.89201375]

>>> o1.symmetry = Oh
>>> o2.symmetry = Oh
>>> print('Oh symmetry misori', np.rad2deg((o1 - o2).angle.data))
>>> print('Oh symmetry angle_with', np.rad2deg(o1.angle_with(o2).data))
Oh symmetry misori [36.81299586]
Oh symmetry angle_with [34.11563997]

from orix.quaternion import Misorientation
>>> misorientation = Misorientation(o1 * ~o2).squeeze()
>>> misorientation.symmetry = (Oh, Oh)
>>> print('Misori', np.rad2deg(misorientation.angle.data))
>>> print('Misori with symmetry reduction', np.rad2deg(misorientation.map_into_symmetry_reduced_zone().angle.data))
Misori [136.89201375]
Misori with symmetry reduction [36.81299586]

In the second run of tests here it looks like angle_with does not compute the same angle as Misorientation. Rotation.angle is written here and looks to be doing a different calculation, but it's a good question as to why the two calculations are producing different values.

[hakonanes]

I can shed some light on the angle_with situation, at the moment the function is written in a way that does not take into account the Orientation.symmetry, perhaps this should be updated!

Symmetry is taken into account via the dot product, so we find the symmetry reduced misorientation angle (disorientation angle)

orix/orix/quaternion/orientation.py

Lines 555 to 568 in 0fa14d1

	def dot(self, other):
	"""Symmetry reduced dot product of orientations in this instance
	to orientations in another instance, returned as
	:class:`~orix.scalar.Scalar`.
	See Also
	--------
	dot_outer
	"""
	symmetry = self.symmetry.outer(other.symmetry).unique()
	misorientation = (~self) * other
	all_dot_products = Rotation(misorientation).dot_outer(symmetry).data
	highest_dot_product = np.max(all_dot_products, axis=-1)
	return Scalar(highest_dot_product)

That said, I wrote this implementation, so I would appreciate a different pair of eyes on it (:

[soupmongoose]

For the final example I made some tests for your findings but was not able to reproduce the angles you quote in the final example:

If I am not mistaken, from_euler should only be given the euler angle in radians not degrees, where you are passing them as degrees, which is why you get different values.

This should work.

o1 = Orientation.from_euler(np.radians([178.5 , 36.97748349 , 63.61592335]),symmetry=symmetry.Oh)
o2 = Orientation.from_euler(np.radians([ 0. , 36.48968436 , 61.39659559]),symmetry=symmetry.Oh)
print('Oh symmetry misori', np.rad2deg((o1 - o2).angle.data))
print('Oh symmetry angle_with', np.rad2deg(o1.angle_with(o2).data))

[harripj]

If I am not mistaken, from_euler should only be given the euler angle in radians not degrees, where you are passing them as degrees, which is why you get different values.

@soupmongoose you are absolutely right, thanks for correcting me!

In any case there is a problem here, the correct value for the symmetry is one thing, but another problem is that the two functions are producing different results.

Looking closer I think I have found a small bug:

o1 = Orientation.from_euler(np.deg2rad([178.5 , 36.97748349 , 63.61592335]))
o2 = Orientation.from_euler(np.deg2rad([ 0. , 36.48968436 , 61.39659559]))

In Orientation.__sub__ the resulting Misorientation is defined such that o2 - o1 returns the transformation from o1 to o2, so the following should return the identity quaternion and does:

# misorientation o2 - o1 is the transformation from o1 -> o2
>>> m = o2 - o1
>>> (m * o1) * ~o2
Orientation (1,) 1
[[ 1. -0. -0.  0.]]

however in Orientation.dot I think that the misorientation is incorrectly defined:

orix/orix/quaternion/orientation.py

Lines 555 to 568 in 0fa14d1

	def dot(self, other):
	"""Symmetry reduced dot product of orientations in this instance
	to orientations in another instance, returned as
	:class:`~orix.scalar.Scalar`.
	See Also
	--------
	dot_outer
	"""
	symmetry = self.symmetry.outer(other.symmetry).unique()
	misorientation = (~self) * other
	all_dot_products = Rotation(misorientation).dot_outer(symmetry).data
	highest_dot_product = np.max(all_dot_products, axis=-1)
	return Scalar(highest_dot_product)

# in Orientation.dot
>>> m = ~o1 * o2
>>> (m * o1) * ~o2
Orientation (1,) 1
[[ 0.8087  0.274  -0.5205 -0.0015]]

which is not identity. This I think is the reason for the two differing results.

In dot I believe disorientation should be defined as misorientation = other * ~self:

>>> m = o2 * ~o1
>>> (m * o1) * ~o2
Orientation (1,) 1
[[ 1. -0. -0.  0.]]

I would like to get your opinions on this @hakonanes @pc494.

NB this is also present in Orientation.dot_outer:

orix/orix/quaternion/orientation.py

Lines 570 to 583 in 0fa14d1

	def dot_outer(self, other):
	"""Symmetry reduced dot product of every orientation in this
	instance to every orientation in another instance, returned as
	:class:`~orix.scalar.Scalar`.
	See Also
	--------
	dot
	"""
	symmetry = self.symmetry.outer(other.symmetry).unique()
	misorientation = (~self).outer(other)
	all_dot_products = Rotation(misorientation).dot_outer(symmetry).data
	highest_dot_product = np.max(all_dot_products, axis=-1)
	return Scalar(highest_dot_product)

I think the line:

orix/orix/quaternion/orientation.py

Line 580 in 0fa14d1

misorientation = (~self).outer(other)

should be:

misorientation = other.outer(~self)

Testing:

>>> m = (~o1).outer(o2)
>>> (m * o1) * ~o2
Orientation (1, 1) 1
[[[ 0.8087  0.274  -0.5205 -0.0015]]]

>>> m = o2.outer(~o1)
>>> (m * o1) * ~o2
Orientation (1, 1) 1
[[[ 1. -0. -0.  0.]]]

[pc494]

Okay, this is slightly off topic, but it's what I spotted in the the skim I was able to do. I think we should deprecated the overridden - op, the difference between two rotations is:

(q1 * ~ q2).angle

writing it other ways is confusing and means there are multiple paths to the same answer.

Will hopefully have time to loop back to this more a more detailed read.

[harripj]

We also should be careful about how we define the misorientation as q2 * ~q1 and q1 * ~q2 are inverse rotations (although the angle will be the same). At the moment o2 - o1 calculates the former and I would favour being consistent with this interpretation.

[hakonanes]

@soupmongoose The misorientation angle between with (o1 - o2).angle is equal to the one calculated with MTEX for all five sets of two Euler angle triplets, even the 39 degree misorientation from (178.5, 36.97748349, 63.61592335) and (0, 36.48968436, 61.39659559). You can continue your work assuming this computation route is reliable, while Orientation.angle_with() needs a patch.

In any case there is a problem here, the correct value for the symmetry is one thing, but another problem is that the two functions are producing different results.

You are right @harripj, thanks a bunch for spotting this bug! After making your changes, the two routes of (o1 - o2).angle and o1.angle_with(o2) return the same results for all five sets of two Euler angle triplets.

We also should be careful about how we define the misorientation as q2 * ~q1 and q1 * ~q2 are inverse rotations (although the angle will be the same). At the moment o2 - o1 calculates the former and I would favour being consistent with this interpretation.

This is another convention we should bring up in #249: if o1 is the reference orientation, i.e. the misorientation is the rotation to bring o1 into o2, o1 should be inverted, as @harripj writes. Checking two text books I have handy (both define orientations from the sample to the crystal reference frame):

• gAB = gB * gA^-1 is written by Stuart Wright (EDAX TSL) in "Texture and Anisotropy" by Kocks, Tomé and Wenk
• M12 = g1^-1 * g2 is written in "Introduction to Texture Analysis" by Engler and Randle

Based on this, I prefer the syntax o2 * ~o1.

---

Scripts comparing orix to MTEX.

orix:

import numpy as np
from orix.quaternion import Orientation, symmetry
eulers = np.deg2rad([
    [[120, 54.26282923, 45], [0, 54.73561032, 45]],
    [[302, -52.35178357, -46.0036662], [0, 53.59948955, 47.5]],
    [[308.5, -51.76150996, -52], [0, 51.76150996, 52]],
    [[328.5, -47.2657904, -67.5], [0, 47.2657904, 67.5]],
    [[178.5, 36.97748349, 63.61592335], [0, 36.48968436, 61.39659559]],
])
for i in range(eulers.shape[0]):
    print(i)
    o1 = Orientation.from_euler(eulers[i][0], symmetry.Oh)
    o2 = Orientation.from_euler(eulers[i][1], symmetry.Oh)

    m12_1 = o1 - o2
    o3_1 = m12_1 * ~o2 * o1
    m12_2 = ~o1 * o2
    o3_2 = m12_2 * ~o2 * o1
    print(m12_1.data.round(4))
    print(m12_2.data.round(4))
    print(o3_1.data.round(4))
    print(o3_2.data.round(4))
    print(np.rad2deg(m12_1.angle.data[0]))
    print(np.rad2deg(o1.angle_with(o2).data[0]))

orix output:

0
[[-1.      0.     -0.0029 -0.0029]]
[[ 0.5    -0.0021 -0.0036  0.866 ]]
[[-5.025e-01  5.000e-04 -5.000e-03  8.646e-01]]
[[ 1.  0. -0.  0.]]
0.47278108586829465
0.47278108586796685
1
[[0.9999 0.0028 0.0082 0.0124]]
[[ 0.0077 -0.4746  0.2721 -0.837 ]]
[[-0.0017  0.4848 -0.2685  0.8324]]
[[ 1. -0.  0.  0.]]
1.7275579455731056
1.727557945573067
2
[[-1.e+00 -6.e-04 -5.e-04 -6.e-04]]
[[ 0.0009 -0.4355  0.2101 -0.8753]]
[[-2.000e-04 -4.361e-01  2.104e-01 -8.750e-01]]
[[ 1.  0. -0.  0.]]
0.11365802066515938
0.11365802066716485
3
[[-1.e+00 -2.e-04 -6.e-04 -6.e-04]]
[[-8.000e-04 -2.705e-01  7.630e-02 -9.597e-01]]
[[ 0.0014 -0.2711  0.0764 -0.9595]]
[[ 1.  0.  0. -0.]]
0.1002851900674705
0.1002851900666658
4
[[-0.9419 -0.1916 -0.1935 -0.1969]]
[[ 0.0024 -0.0116 -0.0041 -0.9999]]
[[ 0.1976 -0.2041  0.1849 -0.9408]]
[[1. 0. 0. 0.]]
39.26841870589627
39.268418705896316

MTEX:

euler = [[[120, 54.26282923, 45], [0, 54.73561032, 45]],...
    [[302, -52.35178357, -46.0036662], [0, 53.59948955, 47.5]],...
    [[308.5, -51.76150996, -52], [0, 51.76150996, 52]],...
    [[328.5, -47.2657904, -67.5], [0, 47.2657904, 67.5]],...
    [[178.5, 36.97748349, 63.61592335], [0, 36.48968436, 61.39659559]]...
] * degree;
euler = reshape(euler, [3, 2, 5]);

for i=1:5
    disp(i - 1)
    cs = crystalSymmetry('m-3m');
    o1 = orientation.byEuler(euler(:, 1, i)', cs);
    o2 = orientation.byEuler(euler(:, 2, i)', cs);

    m12 = inv(o1) * o2;
    o3 = inv(o2) * o1 * m12;
    disp([o3.phi1, o3.Phi, o3.phi2])
    angle(o1, o2) / degree
end

MTEX output:

     0
     0     0     0
ans =
    0.4728

     1
     0     0     0
ans =
    1.7276

     2
     0     0     0
ans =
    0.1137

     3
     0     0     0
ans =
    0.1003

     4
     0     0     0
ans =
   39.2684

[hakonanes]

Prematurely hit "Comment", edited above!

[harripj]

@hakonanes thanks for the detailed analysis and clarification!

Based on this, I prefer the syntax o2 * ~o1.

I agree and in fact I think the other method (Engler and Randle) does not work for us here, which I believe is a passive vs active case?

>>> o1, o2 = Orientation.random((2,))
>>> o2 * ~o1, ~o1 * o2
(Orientation (1,) 1
 [[-0.0474 -0.8901  0.0946 -0.4433]],
 Orientation (1,) 1
 [[-0.0474 -0.3106  0.9068 -0.281 ]])

which are not the same.