get_beam_directions_grid() does not sample fundamental region of hexagonal systems

[emichr]

Describe the bug When using get_beam_directions_grid with diffsims-0.5.dev the entire fundamental region of hexagonal systems is not sampled. With 5 degree sampling resolution, only 162 orientations are sampled, while 386 orientations was sampled with diffsims-0.4.2. It seems that this only happens for hexagonal systems (6/mmm). The missing orientations are quite clear from these plots:

Result from diffsims-0.5.dev:

Result from diffsims-0.4.2:

Also, the sampling grid appears different as well, which is not so clear in the figures above. I am not certain which sampling "pattern" is correct, but the sampling in diffsims 0.5.dev appears more "aesthetic".

To Reproduce

from orix import plot
from orix.quaternion import Orientation, Rotation, symmetry
from orix.vector import Vector3d
from diffsims.generators.rotation_list_generators import get_beam_directions_grid
import matplotlib.pyplot as plt
import numpy as np

resolution = 5.0

systems = {
    'cubic': symmetry.Oh,
    'hexagonal': symmetry.D6h,
    'trigonal': symmetry.D3d,
    'tetragonal': symmetry.D4h,
    'orthorhombic': symmetry.D2h,
    'monoclinic': symmetry.C2h, 
    'triclinic': symmetry.Ci
}

fig = plt.figure(figsize=[4*len(systems), 3])

for i, system in enumerate(systems):
    euler_angles = get_beam_directions_grid(system, resolution, mesh='spherified_cube_edge') #Get a uniform sampling of euler angles
    print(f'Simulated {len(euler_angles)} orientations for {systems[system].name}')

    orientations = Orientation(Rotation.from_euler(np.deg2rad(euler_angles)), symmetry=systems[system])
    ax = fig.add_subplot(1, len(systems), i+1, 
                         projection='ipf',
                         symmetry=systems[system],
                         direction=Vector3d.zvector(),
                         hemisphere='upper')
    ax.set_title(systems[system].name, fontweight='bold')
    ax.scatter(orientations)
plt.tight_layout()

[din14970]

The "gridding" methods available don't work very well in the case of hexagonal crystals. My suspicion is that this is related to rounding issues. I think sampling should be a functionality in orix, not diffsims, so I'm moving it in this PR: https://github.com/pyxem/orix/pull/334 where I also added a gridding method with 6-fold symmetry to avoid these issues.

[hakonanes]

The deleted comment was a duplicate.

Thanks for raising this @emichr, I’ll have a look at this before releasing diffsims.

[emichr]

That sounds great @din14970! I assume it will be a few weeks until that will be ready and available?

By the way @hakonanes, the result is the same with Orix 0.8.2, 0.9.0, and 0.9.0.post0.

[hakonanes]

By the way @hakonanes, the result is the same with Orix 0.8.2, 0.9.0, and 0.9.0.post0.

Thanks for pointing this out @emichr. I initially thought that the cause of the different hexagonal sampling was the change of the meaning of Rotation.from_euler(..., direction=crystal2lab) in orix from 0.8.2 to 0.9.0, but according to your comment, and also my own verification of this, this is not the case.

With diffsims 0.4.2 and orix 0.8.2 or 0.9.0

>>> from diffsims.generators.rotation_list_generators import get_beam_directions_grid
>>> rot = get_beam_directions_grid("hexagonal", 2)
>>> rot.shape
(2255, 3)

With diffsims 0.5.0 and orix 0.8.2 or 0.9.0

>>> from diffsims.generators.rotation_list_generators import get_beam_directions_grid
>>> rot = get_beam_directions_grid("hexagonal", 2)
>>> rot.shape
(915, 3)

I checked the rotations returned from the other crystal systems with the four combinations of versions: taking diffsims 0.4.2 and orix 0.8.2 as the ground truth, the rotations are the same except for the hexagonal system. We should get to the bottom of this before releasing diffsims 0.5.0.

[hakonanes]

Solution: Change hexagonal corner from ~(2, -1, -1, 0)~(1, 0, -1, 0) to (9, 1, -10, 0).

In diffsims 0.4.2 the hexagonal sector corners were [(0, 0, 0, 1), (1, 0, -1, 0), (-1, 2, -1, 0)], set by me, which corresponds to the fundamental sector of D6. Since then, @din14970 updated them [(0, 0, 0, 1), (1, 0, -1, 0), (2, -1, -1, 0)] in https://github.com/pyxem/diffsims/commit/318dfccdf8c4a046c223e6fad635513dd08f0ca3, which was ment to correspond to the fundamental sector of D6h, hence the lower number of vectors. The problem with ~[2, -1, -1, 0]~[1, 0, -1, 0] is that it is a rounded, nicer representation of the actual corner [9, 1, -10, 0], corresponding to the vector with azimuthal angle pi / 6. ~[2, -1, -1, 0]~[1, 0, -1, 0] has an azimuthal angle lower than pi / 6, hence the sampling stops before the corner is reached.

Changing the final corner to [9, 1, -10, 0] gives

>>> import numpy as np
>>> from orix.crystal_map import Phase
>>> from orix.quaternion import Rotation, symmetry
>>> from orix.vector import Miller, Vector3d
>>> from diffsims.generators.rotation_list_generators import get_beam_directions_grid
>>> rot_arr = get_beam_directions_grid("hexagonal", 2)
>>> rot = Rotation.from_euler(np.deg2rad(rot_arr))
>>> rot.size
1050
>>> sector = symmetry.D6h.fundamental_sector
>>> v = rot * Vector3d.zvector()
>>> fig = v.scatter(return_figure=True)
>>> fig.axes[0].plot(sector.edges, zorder=5, color="r")
>>> m = Miller(uvw=sector.vertices.data, phase=Phase(point_group=symmetry.D6h))
>>> m
Miller (3,), point group 6/mmm, uvw
[[1.    0.    0.   ]
 [0.    0.    1.   ]
 [0.866 0.5   0.   ]]
>>> m.coordinate_format = "UVTW"
>>> m
Miller (3,), point group 6/mmm, UVTW
[[ 0.6667 -0.3333 -0.3333  0.    ]
 [ 0.      0.     -0.      1.    ]
 [ 0.4107  0.0447 -0.4553  0.    ]]
>>> m.round()
Miller (3,), point group 6/mmm, UVTW
[[  2.  -1.  -1.   0.]
 [  0.   0.  -0.   1.]
 [  9.   1. -10.   0.]]

[hakonanes]

This line in orix explains why the D6h corners are displayed as they are: axes = m.round(max_index=2).coordinates.astype(int)

[hakonanes]

@emichr, I changed the hexagonal corner in #184: could you redo your checks in the top comment from that branch and confirm that the results are as you expect?

You can directly install my branch with python3 -m pip install -e git+https://github.com/hakonanes/diffsims.git@fix-183#egg=diffsims (as per the Python docs).

[hakonanes]

This is now fixed in #184.