Improving and enhancing the creation of beam grid

[din14970]

---

name: Enhancing beam grid
about: An answer to #129

Release Notes

new features, improvement and bug-fix
Summary:

• changed the way in which get_beam_directions_grid works
  • Now accepts the three arguments: crystal_system, resolution, and mesh
  • Mesh accepts 5 options related to how the sphere is meshed: uv_sphere, normalized_cube, spherified_cube_corner (default), spherified_cube_edge, and icosahedral. The generation of the full sphere mesh is broken up into 3 functions: get_uv_sphere_mesh_vertices, get_cube_mesh_vertices, get_icosahedral_mesh_vertices. The kwargs in get_cube_mesh_vertices allows one to create a normalized grid, or a spherified grid in two ways (detailed explanation in docstrings). The icosahedral mesh was copied and adapted from meshzoo. The full mesh is "cropped" as it was before.
  • It is ensured that for all meshes, resolution means the maximum angular distance to the first nearest neighbor from any grid point.
• tests were updated to reflect the change in get_beam_directions_grid. I also test that the maximum angle between nearest neighbor points is smaller than resolution over the entire grid.
• minor bug fix in UV sphere mesh (original meshing function) to ensure the maximum angular deviation to the first nearest neighbor is smaller or equal than resolution.
• minor bug fix in get_beam_directions_grid: in the "cropping" part I round the dot product to 7 decimals to ensure that numerical errors don't exclude points where the dot product is exactly 0.
• added a few utility functions which are used in the meshing functions:
  • normalize_vectors takes a (N,D) array representing N D-dimensional vectors and returns a (N, D) array of those vectors normalized to unit length.
  • _compose_from_faces utility function adapted from meshzoo to construct the icosahedral grid
  • _get_first_nearest_neighbors from a (N,D) array representing a D-dimensional point cloud of N points, return an (N, D) array representing the first nearest neighbor point of each point
  • _get_angles_between_nn_gridpoints if the 3D point cloud can be represented as unit vectors, they are points on a unit sphere and we can calculate the angles between these vectors. The first nearest neighbor of each gridpoint is the vector with the smallest relative angle. This function returns the angles between each grid point and its nearest neighbor. This function could be useful for getting an angular spread distribution for different meshes.
  • _get_max_grid_angle get the maximum angle between any point and its nearest neighbor in a sphere grid
  • _get_min_grid_angle get the minimum angle between any point and its nearest neighbor in a sphere grid
  • _vector_grid_to_euler convert a (N, 3) array representing N 3D vectors to orientations with two euler angles in the bunge convention: (0, Phi, phi2). The orientation in this sense means the rotation to turn [001] to be parallel to the vector, with the first euler angle constrained (or the other way around, I always get confused). I use this function to verify my output with MTEX.

What does this PR do? Please describe and/or link to an open issue.
Previously, the user had limited control in creating a grid of directions on the fundamental triangle. Only a UV-sphere type grid was possible, which significantly oversamples near the (100) poles in the azimuthal direction. The two options (equal angle versus equal area) were non-sensical. A perfect grid on the sphere is not possible, but more fair options exist, each with their own drawbacks. To give the user maximum control, I added 4 new methods of gridding the sphere (described above). In each case, I ensure that the resolution argument is the maximum angle from any point on the grid to its nearest neighbor.
See more detailed descriptions in #129 and #128, as well as pyxem/orix#125

Describe alternatives you've considered
All alternatives I could think of I implemented. One could still go for more exotic sphere meshes but as the resolution is refined, the returns will probably be minimal.

Are there any known issues? Do you need help?
I tried to iron out the issues I could find. One could still verify the corners of the fundamental zone as there were some discussions about this in #128

Work in progress?
If you know there is more to do, make a checklist here:

• I now implemented everything in the file rotation_list_generator for simplicity. However, I can imagine that some helper functions do not belong there and may better fit elsewhere. This I leave up to you.
• Optional: Testing could always be improved, utility functions are tested implicitly. One should test for the "correctness" of each grid. Now I only test for maximum angle.
• Optional: The icosahedral grid is slow because it is generated iteratively, i.e. a number of starting subdivisions is chosen, the maximum angle on the grid is calculated, if it is larger than resolution, the number of subdivisions is increased. If a maximum angle could be known a priory, the generation of the grid would be much faster.

[pc494]

This looks pretty good to me at a first glance; it'll probably take two go rounds with the review but I'll try and get the first one in at some point this week.

[pc494]

This is a good PR, and most of these requests are for/about readability or clarity. I've found in the long term such changes pay off.

[din14970]

still found some bugs, will fix now

[pc494]

This looks good to me, I will make a couple of packaging additions and then merge.

[pc494]

Sorry @din14970, I have another question. I was looking at beam_directions_grid_to_euler and hoping to merge it into get_rotation_from_z_to_direction (from zap_map_generator.py) but I think there might be a disagreement about euler angle conventions.

diffsims would treat bunge euler angles from left to right and so I would expect phi2 (the final angle) to be the zero and not phi1?

[din14970]

@pc494 This also always messes me up and has to do with active vs. passive rotations and also how the euler angles are defined. Actually thanks for making me have a second look at this because it appears what I have is still not quite right.

Just to clarify my solution:

what I looked at was what I got out of the MTEX plots and what information I can derive from the ASTAR software. So I verified for the cubic case that the rotations I export with my functions do give me back the right IPF plots in MTEX. In the cubic case we want to populate the triangle [0, 0, 1] - [1, 0, 1] - [1, 1, 1]. I verified that this works properly when generating the grid. I then convert to Euler angles with beam_directions_grid_to_euler, export, and import to MTEX. I would then expect the full IPF to have this triangle populated with points, because the IPF-Z should tell me which crystal vectors to expect parallel to the world Z-axis in my collection of orientations. The IPF images look like:

So while 001 and 111 are ok, I have the wrong 011, which I obviously overlooked because the shape looked ok. Changing phi2 to np.pi/2 - phi2 produced the expected triangle:

Maybe the right way to do it requires a bit more thought and proper testing.

Regardless, flipping phi1 and phi2 results in the following IPF, which is clearly incorrect:

A final bit of evidence is a screenshot from the ASTAR indexing software:

Clearly it is the phi1 angle that is optimized on indexing through the polar transform - the graph then shows the correlation as a function of phi1. Also the template bank has all templates with phi1 = 0.

TL:DR:
I'm quite confident it is phi1 that should be zero, however I think I'm off by a shift of pi/2 on phi2 for some reason.

[din14970]

As a small note, my function adjusted with the pi/2 term is also able to reproduce the ASTAR angles in the last image from the previous comment:

I think the error is from rounding errors in trying to make the miller indices integers.

[pc494]

I think I was being a little slow, because there is a spare degree of freedom we may well be using the same convention and getting "different" answers that both produce the same crystal normal.

If you don't mind, I think the next steps ought to be

1. check (in python) that the beam directions (cartesian vectors) fill the correct streogram after gridding (fairly sure they do)
2. fix get_rotation_from_z_to_direction so that you can take a grid (of vectors) and export them as a grid of euler angles to import correctly to MTEX

The reason for doing it like this is that writing code that makes use of euler angles has always seemed very risky. Especially when we are working in a fairly simple geometry ([1] captures my feelings nicely).

[1] - https://github.com/moble/quaternion/wiki/Euler-angles-are-horrible#opening-issues-and-pull-requests

[din14970]

Are you proposing I do this for this PR?

check (in python) that the beam directions (cartesian vectors) fill the correct streogram after gridding (fairly sure they do)

you mean with some kind of test on the correctness of get_beam_directions_grid?

fix get_rotation_from_z_to_direction so that you can take a grid (of vectors) and export them as a grid of euler angles to import correctly to MTEX

Yes obviously if there's already a function for doing this then my beam_directions_grid_to_euler is unnecessary. But correct me if I'm wrong, but it would seem that get_rotation_from_z_to_direction gives the "active" rotation from z onto the crystal direction, whereas I've read that crystal orientations are always defined in the "passive" way. So should the function not be get_rotation_from_direction_to_z? I'll check what I get in MTEX.

An additional comment: I noticed that you need a crystal structure in get_rotation_from_z_to_direction which is obviously the case when "direction" represents a crystal direction [uvw] or (hkl). However, in the case of the grid, one should not interpret these directly as crystal vectors. Variation among lattice parameters will mean that you significantly distort the grid that way. For example, if you have a tetragonal material with say a=1 and c=10, then the number of grid points between [001] and [101] will be the same as the number of gridpoints between [100] and [101], but the angular deviation between the latter will be much greater. So the grid of vectors should be converted to euler angles without taking into consideration the underlying lattice.

Perhaps this calls for having a separate grid conversion function?

[pc494]

Are you proposing I do this for this PR?

I can do it. I will just plot the streogram directly using the formulae from wikipedia.

fix get_rotation_from_z_to_direction so that you can take a grid (of vectors) and export them as a grid of euler angles to import correctly to MTEX:

Yes obviously if there's already a function for doing this then my beam_directions_grid_to_euler is unnecessary.

Change of heart. I think I got up to speed on what you're trying to do (bigger picture). I had always expected to do template matching in two phases. Radial. Get results. In plane. But you are planning on doing it in a single swoop, which means you're euler angles need to do "good" things (like be 0 at an inplane direction or whatever) - which your function does (or would do if it's broken?!) and the current one doesn't.

But correct me if I'm wrong, but it would seem that get_rotation_from_z_to_direction gives the "active" rotation from z onto the crystal direction, whereas I've read that crystal orientations are always defined in the "passive" way. So should the function not be get_rotation_from_direction_to_z? I'll check what I get in MTEX.

I think there are two separate things. "passive" and "active" rotations tell us whether our axes travel with the rotation of the object (all our code is "passive", the axes move with the system), see for example section 2.3 of [1]. What you seem to be talking about is if we define rotation as "crystal2lab" or "lab2crystal". This has been brought to my attention before, very broadly SPED leans one way and EBSD the other & the two are related by an inversion.

An additional comment: I noticed that you need a crystal structure in get_rotation_from_z_to_direction which is obviously the case when "direction" represents a crystal direction [uvw] or (hkl). However, in the case of the grid, one should not interpret these directly as crystal vectors. Variation among lattice parameters will mean that you significantly distort the grid that way. For example, if you have a tetragonal material with say a=1 and c=10, then the number of grid points between [001] and [101] will be the same as the number of gridpoints between [100] and [101], but the angular deviation between the latter will be much greater. So the grid of vectors should be converted to euler angles without taking into consideration the underlying lattice.

Yeah, I had noticed this too, I had been planning on adding an extra keyword argument "cubic" to treat the input as cartesian (which they are here). I think I will still do that but...

Perhaps this calls for having a separate grid conversion function?

...I also agree with this.

[1] - https://iopscience-iop-org.ezp.lib.cam.ac.uk/article/10.1088/0965-0393/23/8/083501/meta
[1b] if the links doesn't work, it's meant to point to "Consistent representations of and conversions between 3D rotations" by Rowenhorst et al. (2015)

[pc494]

Checked the triangle, looks fine to me, merging.

g = get_beam_directions_grid("cubic",2)
from matplotlib import pyplot as plt

X = np.divide(g[:,0],1-g[:,2])
Y = np.divide(g[:,1],1-g[:,2])

plt.scatter(X,Y)

conversion formulae from:
https://en.wikipedia.org/wiki/Stereographic_projection

[pc494]

This is an excellent addition, many thanks @din14970