Mean of orientations, with symmetry

[DorianDepriester]

Hello there!
It seems that the mean() method for orientations actually doesn't take the orientations' symmetry into account. For instance, let's consider the m-3m symmetry and a set of 3 close orientations:

from orix.quaternion import symmetry, Orientation

s = symmetry.Oh # Cubic symmetry
o = Orientation.from_euler([[0,0,0], [1,89,0], [1,1,89]],symmetry=s,degrees=True)
print(o.in_euler_fundamental_region()*180/np.pi) # Just to be sure what we're talking about

# Compute and print the mean orientation
obar = o.mean()
print(obar.in_euler_fundamental_region()*180/np.pi)

Thus, the output is:

[[ 0. 0. 0.]
[ 1. 89. 0.]
[ 1. 1. 89.]]
[[15.79096054 30.29090638 15.77845771]]

The result is quite strange. Indeed, the same computation on MTEX

CS=crystalSymmetry('m-3m')
mean(orientation.byEuler([0,1,1]*degree, [0,89,1]*degree, [0,0,89]*degree,CS))

gives

ans = orientation (m-3m -> xyz)
 
  Bunge Euler angles in degree
     phi1        Phi       phi2       Inv.
  90.8333 0.00290889      269.5          0

Still, the same with CS=crystalSymmetry('1') leads to:

Bunge Euler angles in degree
phi1 Phi phi2 Inv.
15.791 30.2909 15.7785 0

Am I doing something wrong? I have looked into MTEX's code, and I understand that the mean() method actually relies on the project2FundamentalRegion function. I have tried to implement it in orix, but without any success.

Thank you in advance.

[hakonanes]

Hi Dorian,

you have indeed encountered a somewhat surprising result, from a user standpoint. Orientation.mean() does not consider symmetry, as you say. I won't quite call it a bug since the method description states that it computes the quaternion mean. But, it is something we should support, and can do quite easily.

Do you get the desired results with this approach?

>>> from orix.quaternion import Orientation, symmetry
>>> s = symmetry.Oh
>>> o = Orientation.from_euler([[0, 0, 0], [1, 89, 0], [1, 1, 89]], symmetry=s, degrees=True)
>>> o2 = o.map_into_symmetry_reduced_zone()
>>> obar2 = o2.mean()
>>> obar2.to_euler(degrees=True)
array([[179.9955333 ,   0.32750355, 180.00737556]])

The above is basically the approach we should use in Orientation.mean().

[DorianDepriester]

Thank you @hakonanes for this hack. Seems good!

[hakonanes]

Glad it solves your immediate problem! Please report back here if it does not or you thought of something else related to this issue...

I suggest we leave this open until we make Orientation.mean() account for symmetry.

[DorianDepriester]

Hello @hakonanes,
I have noticed something weird in the code above. For instance, if you try this:

o =  Orientation.from_euler([[80, 42, 10]], degrees=True, symmetry=s)
o2 = o.map_into_symmetry_reduced_zone()
obar2 = o2.mean()
print(obar2.to_euler(degrees=True))

os = Orientation((o.data[0],obar2.data[0]),symmetry=o.symmetry)
print(os.get_distance_matrix(degrees=True))

you get:

[[350.  42.  10.]]
[[ 0.         54.16431688]
 [54.16431688  0.        ]]

As you can see here, the mean of a single orientation is not equal to this orientation. This is evidenced by the non-zero misorientation angle between o and its own mean. Am I doing or understanding something wrong?

Rgds

[hakonanes]

Very good point. I've made a change to map_into_symmetry_reduced_zone() in #442 (also exchanging the name for reduce(), with the former method to be removed two releases from now in v0.13) which solves your obvious issue here:

from orix.quaternion import Orientation, symmetry


s = symmetry.Oh
o =  Orientation.from_euler([80, 42, 10], degrees=True, symmetry=s)
o2 = o.reduce()
obar2 = o2.mean()
print(obar2.to_euler(degrees=True))
# [[ 80.  42. 280.]]

os = Orientation((o.data[0], obar2.data[0]), symmetry=o.symmetry)
print(os.get_distance_matrix(degrees=True))
# [[0. 0.]
# [0. 0.]]

I would be very grateful if you could check out the branch in that PR and try out reduce() in your current use of orix, to see if reduction to the fundamental zone gives more intuitive results than before.

[DorianDepriester]

Hello @hakonanes ,
Thank you for this quick fix. I will have a look on the PR as soon as I enough free time (not today, neither the next 14 days, to be honest…).

[hakonanes]

No problem, I don't think we will release this change within two weeks anyway.

[DorianDepriester]

Hello @hakonanes ,
With the aid of your comment in #442 (comment), I can easily reproduce the results you get. In addition, my former implementation of the weighted mean orientation #435 (comment) seems to work as well. Here is a simple test:

from orix.quaternion import Quaternion,Orientation, symmetry
import numpy as np

def mean_orientation(o, weights=None):
    if weights is None:
        weights=np.ones((1,o.shape[0]))
    else:
        # Force weights to be a 2d array
        weights=np.array(weights, ndmin=2)
    o2 = o.reduce()
    q = o2.data
    qq = np.einsum('pi,ij,ik->pjk', weights, q, q)
    w, v = np.linalg.eig(qq)
    w_max = np.argmax(w, axis=1)
    q_mean= v[np.arange(weights.shape[0]), :, w_max]
    return Orientation(Quaternion(q_mean), o.symmetry)


ori = Orientation.from_axes_angles(
    [[1, 0, 0], [1, 0, 0], [0, 1, 0], [0, 1, 0], [0, 0, 1], [0, 0, 1]],
    [30, -30, 30, -30, 30, -30],
    symmetry=symmetry.Oh,
    degrees=True
)
print(ori.to_euler(degrees=True))
weights=[[1, 1, 1, 1, 1, 1],[1000, 1, 1, 1, 1, 1],[1, 1000, 1, 1, 1, 1],[1,1,1000,1,1,1]]

# Compute weighted mean orientation, wrt. the weighting matrix
ori_w = mean_orientation(ori, weights=weights)

# Print each weighted mean and its misorientation wrt. the initial orientations
for i,w in enumerate(weights):
    print(' ')
    print('Weights: {}'.format(w))
    print('Mean Euler angles: {}'.format(ori_w[i].to_euler(degrees=True)))
    os = Orientation(np.vstack((ori.data,ori_w.data[i])),symmetry=ori.symmetry)
    d=os.get_distance_matrix(degrees=True)
    print('Misorientation angles from weighted mean:')
    print(d[-1,:])

[[180.  30. 180.]
 [  0.  30.   0.]
 [270.  30.  90.]
 [ 90.  30. 270.]
 [330.   0.   0.]
 [ 30.   0.   0.]]

Weights: [1, 1, 1, 1, 1, 1]
Mean Euler angles: [[0. 0. 0.]]
Misorientation angles from weighted mean:
[30. 30. 30. 30. 30. 30.  0.]

Weights: [1000, 1, 1, 1, 1, 1]
Mean Euler angles: [[180.          29.84404743 180.        ]]
Misorientation angles from weighted mean:
[ 0.15595257 30.15595257 42.07295755 42.07295755 42.07295755 42.07295755
  0.        ]

Weights: [1, 1000, 1, 1, 1, 1]
Mean Euler angles: [[ 0.         29.84404743  0.        ]]
Misorientation angles from weighted mean:
[30.15595257  0.15595257 42.07295755 42.07295755 42.07295755 42.07295755
  0.        ]

Weights: [1, 1, 1000, 1, 1, 1]
Mean Euler angles: [[270.          29.84404743  90.        ]]
Misorientation angles from weighted mean:
[42.07295755 42.07295755  0.15595257 30.15595257 42.07295755 42.07295755
  0.        ]

Great job!

[hakonanes]

Good! Thank you for doing the tests. Maintainer resources are in short supply, so we just have to wait a bit for #442 to be reviewed.

[argerlt]

Just to add to this:

mtex doesn't just put points into the fundamental zone to find the mean. It starts with that, but then it finds the approximate mean, calculates a new fz with it's center aligned with the roughly calculated orientation mean,, re-pushes everything to the new fz, then does a final mean calculation with the new set of points.

It does this because otherwise, if the cluster happens to sit close to a boundary, the Euclidean mean is incorrectly calculated because some points are on the opposite side of a non-mirror symmetry operation, and the answer gets skewed.

like @hakonanes said, it's technically not an error since orix only claims to calculate the quaternion mean, but some day we should add a symmetry-aware mean.

As a fun aside, the mtex solution is better but still imperfect, and a proper solution is one of those astoundingly difficult problems to solve that has led to some truly wild solutions. My favorite example is the MTEX embeddings method based on this paper, which is able to find accurate means, but only by embedding orientations into symmetry specific tensors to get 10+dimensional vector representations.

Non-Euclidean math is wild.

[hakonanes]

Yes, we can implement the second pass for calculating the mean. Shouldn't be a problem.