Option to use euler fundamental region for mean orientation

anderscmathisen commented 1 year ago

Description of the change

Added override to the mean() function in the Orientation class with the option to calculate the mean in the euler fundamental region. Also keeps the symmetry of the Orientation for the returned mean. This way, when plotting in an IPF, I think the mean Orientation makes more sense.

Progress of the PR

[x] Docstrings for all functions
[x] Unit tests with pytest for all lines
[x] Clean code style by running black via pre-commit

Minimal example of the bug fix or new feature

import orix.plot
from orix.quaternion import Orientation
from orix.quaternion.symmetry import D6h
from orix.vector import Vector3d
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec

oris = Orientation(
[[ 0.215 ,  0.9696, -0.1001, -0.0602],
 [ 0.2132 , 0.97  , -0.0977, -0.0643],
 [ 0.2154 , 0.9693 ,-0.098 , -0.067 ],
 [ 0.2157 , 0.9692, -0.0989, -0.0659],
 [ 0.2141 , 0.9693 ,-0.1012, -0.0668],
 [ 0.2132 , 0.9696, -0.1001, -0.0659],
 [ 0.2083 , 0.9715 ,-0.0952, -0.0605],
 [ 0.1558 , 0.7913 ,-0.569 , -0.1608],
 [ 0.1495 , 0.79   ,-0.5704, -0.1681],
 [ 0.1537 , 0.7908, -0.5686, -0.1663],
 [ 0.1552 , 0.7918 ,-0.5677, -0.1634],
 [ 0.1518 , 0.792 , -0.5679, -0.165 ],
 [ 0.0484 , 0.3993, -0.8886, -0.2203],
 [ 0.0459 , 0.4007, -0.888,  -0.2209]], symmetry = D6h)

fig = plt.figure(figsize=(9,5))

gs = GridSpec(2,3)

ax_x_not_fundamental = fig.add_subplot(gs[0,0], projection = "ipf", direction = Vector3d.xvector(), symmetry = D6h)
ax_y_not_fundamental = fig.add_subplot(gs[0,1], projection = "ipf", direction = Vector3d.yvector(), symmetry = D6h)
ax_z_not_fundamental = fig.add_subplot(gs[0,2], projection = "ipf", direction = Vector3d.zvector(), symmetry = D6h)

ax_x_fundamental = fig.add_subplot(gs[1,0], projection = "ipf", direction = Vector3d.xvector(), symmetry = D6h)
ax_y_fundamental = fig.add_subplot(gs[1,1], projection = "ipf", direction = Vector3d.yvector(), symmetry = D6h)
ax_z_fundamental = fig.add_subplot(gs[1,2], projection = "ipf", direction = Vector3d.zvector(), symmetry = D6h)

ax_x_not_fundamental.set_title("X")
ax_y_not_fundamental.set_title("Not in fundamental region\nY")
ax_z_not_fundamental.set_title("Z")

ax_x_fundamental.set_title("X")
ax_y_fundamental.set_title("In fundamental region\nY")
ax_z_fundamental.set_title("Z")

mean_ori_not_in_fundamental = oris.mean()
mean_ori_in_fundamental = oris.mean(mean_in_euler_fundamental_region=True)

ax_x_not_fundamental.scatter(oris, label ="Data")
ax_y_not_fundamental.scatter(oris)
ax_z_not_fundamental.scatter(oris)
ax_x_not_fundamental.scatter(mean_ori_not_in_fundamental, c = "r", label = "Mean")
ax_y_not_fundamental.scatter(mean_ori_not_in_fundamental, c = "r")
ax_z_not_fundamental.scatter(mean_ori_not_in_fundamental, c = "r")

ax_x_fundamental.scatter(oris)
ax_y_fundamental.scatter(oris)
ax_z_fundamental.scatter(oris)
ax_x_fundamental.scatter(mean_ori_in_fundamental, c = "r")
ax_y_fundamental.scatter(mean_ori_in_fundamental, c = "r")
ax_z_fundamental.scatter(mean_ori_in_fundamental, c = "r")

fig.tight_layout()
handles, labels = ax_x_not_fundamental.get_legend_handles_labels()
handle_list, label_list = [], []
for handle, label in zip(handles, labels):
    if label in ["sa_circle", "sa_sector"]:
        continue
    if label not in label_list:
        handle_list.append(handle)
        label_list.append(label)
fig.legend(handle_list, label_list, loc = "upper left", fontsize = 12)

Orientation

For reviewers

[ ] The PR title is short, concise, and will make sense 1 year later.
[ ] New functions are imported in corresponding __init__.py.
[ ] New features, API changes, and deprecations are mentioned in the unreleased section in CHANGELOG.rst.
[ ] Contributor(s) are listed correctly in __credits__ in orix/__init__.py and in .zenodo.json.

hakonanes commented 1 year ago

Thanks for this, @anderscmathisen.

Did you consider using Orientation.map_into_symmetry_reduced_zone().mean() instead? I would expect this to give correct results. But then I checked, and it does not... I knew there were some inconsistencies between this function and the other projection tools in orix. This has been bothering me for a while, and I believe we should double check this implementation before continuing with this PR.

I actually have a rewritten version of EMsoft's (Singh and De Graef's) reduction to the fundamental zone, and it produces the same results as you get by projecting to the Euler fundamental region. Interestingly, it reduces all your orientations above to a tight cluster in the axis-angle space of 622. Is this as expected?

For completeness' sake, this change also relates to #434 in that we should overwrite more methods in the Orientation class to better handle properties (like symmetry).

anderscmathisen commented 1 year ago

Thanks for the swift response.

I indeed initially tried to use map_into_symmetry_reduced_zone(), but when it did not provide the desired results, I choose the euler fundamental zone instead.

I agree that the map_into_symmetry_reduced_zone() would be more natural to use here. I dont yet really understand how that function is supposed to work, and thus have not figured out why it does not. When plotting oris from the code above in 3D (oris.map_into_symmetry_reduced_zone().scatter()) the orientations appear to be left untouched e.g. same as oris.scatter(). It appears to me that Misorientation.map_into_symmetry_reduced_zone() is coded to work best with Misorientations and perhaps not the child Orientation? A possible (quick) fix is to overwrite that function in the Orientation class using in_euler_fundamental_region()? But map_into_symmetry_reduced_zone() is perhaps supposed to be fundamentally different from in_euler_fundamental_region()? I expected them to be the same.

I actually have a rewritten version of EMsoft's (Singh and De Graef's) reduction to the fundamental zone, and it produces the same results as you get by projecting to the Euler fundamental region. Interestingly, it reduces all your orientations above to a tight cluster in the axis-angle space of 622. Is this as expected?

Not really no.

Edit: My implementation here does not fix the example problem in #434, so I guess something is wrong and map_into_symmetry_reduces_zone()should be used. However this would not work for my example here...

hakonanes commented 1 year ago

Not really no.

OK, I checked my results against EMsoft's using their command line program EMConvertOrientations, converting quaternions to quaternions using point group 6/mmm and reducing them to the fundamental zone. EMsoft gives three clusters, as the current implementation in orix (although they are different...). My implementation misses something, which I'll fix. We can hopefully then update the current implementation of map_into_symmetry_reduced_zone(). This should then fix both your issue here and the one in #434.

hakonanes commented 1 year ago

I've changed reduction to the Rodrigues fundamental zone in #442. This change should fix your issue here, i.e. if you pull that branch and replace

mean_ori_in_fundamental = oris.mean(mean_in_euler_fundamental_region=True)

with

mean_ori_in_fundamental = oris.reduce().mean()

you should get what you want.

I would be very grateful if you could actually check out the change on your example here and also any other problems you might have where reducing to the fundamental zone would be desirable.

anderscmathisen commented 1 year ago

Hi,

I have been testing with #442 and indeed that seems to solve all the issues with mean orientation :)

As we discussed offline @hakonanes, we should rebase this PR to #442 and replace my previously committed code with something along the lines of:

def mean(self, reduce=True):
    """..."""
    if reduce:
        quat = Quaternion(self.reduce())
    else:
        quat = Quaternion(self)
    return Orientation(quat.mean(), symmetry = self.symmetry)

It is perhaps best to wait for #442 to merge into develop before rebasing and continuing with this PR?

hakonanes commented 1 year ago

Great. Yes, since we might make some changes to #442 before merging.

hakonanes commented 1 year ago

BTW, since Misorientation.reduce() will be updated to align with expectations in #442 as well, I suggest to upstream your change here to Misorientation.mean():

def mean(self, reduce=True):
    """..."""
    if reduce:
        quat = Quaternion(self.reduce())
    else:
        quat = Quaternion(self)
    return self.__class__(quat.mean(), symmetry=self.symmetry)

pyxem / orix