Fixing get_grid_beam_directions

[pc494]

---

name: Fixing get grid directions
about: In light of pyxem/orix#125 & #103 it seems get_grid_beam_directions was producing bad answers, this PR attempts to fix that.

---

Release Notes
See changelog.

What does this PR do? Please describe and/or link to an open issue.
The routine used has been changed. It is now:

1. Cover the whole sphere with directions (angle or area supported)
2. Remove the duplicates that occurs at theta=0 and theta=180
3. Keep only those direction that lie in the relevant triangle via the following:

Take the 3 corners of the triangle.

a) Each pair of corners forms 2 vectors (call them A,B) with the origin that define a plane, find normal of said plane: N
b) You only want points that lie on the side of the third corner (C), which are vectors with dot products of the same sign as (c dot N)
c) Keep those points, remove all the other ones
d) Repeat with A,C and B,A

Return the relevant cartesian coordinates

Help wanted
Testing ideas or thoughts are very welcome, cc: @hakonanes and @din14970

[hakonanes]

Pulled this PR and ran @din14970's code bit from pyxem/orix#125 (comment) with a step of 1 degree and both angle and area. Equal angle/area produced 2181/2733 rotations. Relevant MTEX lines:

cs = crystalSymmetry.load(fullfile(cifdir, '1534888.cif'));  % m-3m

% Read orientations
o = orientation.load(file, 'columnNames', {'phi1', 'Phi', 'phi2'});
o2 = orientation.load(file, 'columnNames', {'phi1', 'Phi', 'phi2'}, cs);

% Plot inverse pole figures
figure
plotIPDF(o2, [xvector, yvector, zvector], 'all')

% Plot full pole figures
figure
plotIPDF(o2, [xvector, yvector, zvector], 'upper', 'all', 'complete')
hold on
plotIPDF(o, [xvector, yvector, zvector], 'upper', 'all', 'MarkerFaceColor', 'r')

Equal angle:

Equal area:

Equal angle:

Equal area:

At a glance this looks good to me, although I need to look a bit more at your code to understand the plots better.

[hakonanes]

@pc494 which source have you used for the triangle corners? Using 11-20 for hexagonal doesn't cover the full pole figure, while -12-10, which MTEX uses, does. Corners 11-20/-12-10 produced, for equal area and 1 degree step, 1692/5697 orientations. Using same MTEX code for IPFs as above, only with crystalSymmetry('6'), produces:

11-20 corner in diffsims (x upper/lower, y upper/lower, z upper/lower):

-12-10 corner in diffsims:

[pc494]

The corners came from an old bit of Simon's code, I struggled to find a more comprehensive source. Taking the ones from MTEX seems like a good way forward.

[hakonanes]

Okay, if it's okay with you @pc494 I'll commit to this PR with the corners MTEX uses when plotting orientations generated by diffsims with the seven different crystal systems. The MTEX crystal symmetries are then created using crystalSymmetry('cubic') (m-3m), crystalSymmetry('hexagonal') (6/mmm) etc.

[pc494]

Yep, although the code currently skips triclinic to just return everything, probably worth making that work in the same way as everything else.

[din14970]

@pc494 I think the cropping based on corners is a quite elegant solution.

However, I'm not completely on board with the grids based on polar coordinates - the points near 001 are very close together in the azimuthal direction, no matter if area or angle is chosen. Would it be possible to also add my "gridding" method based on gridding a cube surrounding the sphere? I find this:

to be fairer looking than the grids in the images from @hakonanes especially near 001.

In addition, I am now even more certain that this is also how Astar does it, and it would be nice at least to have this option available for direct comparison. For the cube case specificaly, to cover only the fundamental triangle I did:

def get_zone_vectors_cube(steps=50):
    k = np.arange(0, steps+1)
    theta = np.arctan(np.sqrt(2))/steps
    i = np.tan(k*theta)/np.sqrt(2)
    x, y = np.meshgrid(i, i)
    x, y = x.ravel(), y.ravel()
    yy = y[y<=x]
    xx = x[y<=x]
    z = np.ones(xx.shape[0])
    all_vecs = np.vstack([xx, yy, z]).T
    return all_vecs

This gave me 1326 points, exactly the same number as the number of templates in my ASTAR bank file:

Note also the correspondence with step count. Instead of using steps as an argument one could change this to a "maximum_angle" argument which would be the angle between 001 and the first grid point towards the 111 direction.

I wouldn't throw away the grid based on spherical coordinates though because maybe it might be useful in other circumstances. Would it be best to add my gridding method as an option in the get_beam_directions_grid function with an additional argument, or add it as a separate function?

[pc494]

@din14970 could you just clarify the comment "the points near 001 are very close together in the azimuthal direction, no matter if area or angle is chosen." for me?

My intention had been that in the "angle" case the directions are all "resolution" away from their 2/4 nearest neighbours?

[din14970]

@pc494 sorry about that, it gets very confusing discussing vectors and orientations with mainly text. A virtual blackboard would be convenient.

I mean that since the grid is based on polar coordinates, the number of points on each "lattitude" line is the same. So near the the [001] pole, the spacing between the points on the same lattitude line is tiny and it appears to be nearly a continuous line

Whereas more towards the "equator" the spacing becomes more apparent

The "resolution" is only valid on the equator; the angular spacing on the same lattitude line is much much smaller. Only the spacing between the lattitude lines remains valid if you choose equal angle.

Bottom line, you are not very fairly sampling the surface of the sphere with this method in my opinion.

[pc494]

I think this might be an actual science disagreement then.

My case would be that the distance (using a globe analogy) that we care about is the angular one and not the "walking on the surface of a sphere one". A tilt of 1 degree in a SED experiment is always 1 degree, regardless of how "far" a surface normal has travelled. I want my grid to mean that if I have 1 degree resolution about an 001 zone axis, I also have 1 degree around the 111 zone axis.

Given that it's science here, I'll get thinking on how best fit what you've written in; does it transfer to other crystal systems cleanly?

EDIT (addition): To illustrate how the two different aims are incompatible.
The gridding you've provided has n points along each of the edges of the stereographic triangle. But the 001-011 edge corresponds to a rotation of 45°, while the 001-111 edge is 54.7°, resulting in what I would describe as "asymmetric (angular) resolution". This effect is most pronounced at the 111 pole, as in one direction you get 54.7/n and in the other you get 35.2/n

EDIT 2 (further remark)
I think this is a long way of saying, "feature not a bug" wrt angle and area, Adding your method is a good idea, can you raise an "enhancement issue" for it, as there are some compatibility things its worth being mindful of. Meanwhile I think I'll merge this once @hakonanes has been able to get his bit in.

[din14970]

@pc494 Maybe I'm totally misunderstanding your point, but the way I read it I think you are glossing over something. Apologies if I misread, if I'm totally off base, or if I'm pedantic. Sorry also for the long comment, I don't know how to explain what I mean more efficiently. Maybe the issue is that you may have a different use case in mind for this type of function?

For me this type of function would serve to generate a list of vectors representing beam directions, from which I want to generate templates. Therefore I want my vectors to sample the fundamental zone as best and fairly as possible. As best as possible means indeed that the average angle between nearest vectors is the same for all vectors. I allude to "distance" and "far" in my previous comment, but this is equivalent to angle because arc length between two points on a sphere surface is proportional to the angle between the lines connecting those points to the center of the sphere.

Unfortunately in the general case it is impossible to completely fairly sample a sphere surface with a regular grid, because a sphere is topologically different from a plane - you can not draw points on a sphere so that the arc length between all points is the same unless those points represent the vertices of a platonic solid. So the question becomes how can one do this to the best approximation.

And so, I agree with you that my solution is suboptimal, as indeed I have 50 grid points on all sides of the stereographic triangle, while the arc-length between each corner is not the same. What I am saying is that a grid based on spherical coordinates is even less fair, no matter if you choose angle or area. You state:

I want my grid to mean that if I have 1 degree resolution about an 001 zone axis, I also have 1 degree around the 111 zone axis.

This is not achieved with a spherical coordinate grid. Suppose I choose equal angle and a 1 degree resolution. Then for all the points on the equator of the sphere, your grid will be completely fair: the angular spacing between 4 nearest neighboring grid points is always 1 degree. But this is true only at the equator; if you move further up the sphere to [001] or down to [00-1], the angular spacing between elevation lines remains 1 degree, but the points on that elevation line get closer together. Only the projected vectors onto the x-y plane are still spaced 1 degree from eachother - this is delta psi; the actual angular difference between the points gets smaller. You can verify this with a dot product calculation - the points near [111] or [011] are much more closely spaced in angle than the points on the equator. The good thing is your 1 degree distance is "worst case", but you are significantly and in my view unnecessarily oversampling near the poles of the sphere. So just like my case, you have "assymetric" resolution.

You try to mitigate this with the area option and you refer to a wolfram page in the source code. But this is only valid for randomly sampling the sphere surface - the point is if you parametrize the sphere with spherical coordinates and you pick points by drawing psi and theta from a uniform distribution you will oversample the poles for exactly the same reason as I point out earlier: points get closer on higher elevation lines. So the solution for this random sampling is to correct by surface area. But if you do this transformation on your grid, all you achieve is having fewer elevation lines near the [001] pole. So the distance between the 4 nearest neighbors of each grid point will be even more unfair - in the azimuthal direction the distance will be always <1 degree depending on how high up you are, but in the elevation direction it will be >1 degree.

Probably the fairest "grid" is a randomly generated one, which is eroded such that the average and maximum angular distances between points meet certain criteria. But then you don't really have much control or reproducibility.

While my grid is not completely fair, as a completely fair grid is impossible, I believe it is more fair in terms of nearest neighbor angular spacing than the spherical coordinate grid. And as you increase the step size, the unfairness in sampling should become less significant. It is of course best suited for a cube as there will always be gridpoints on the 110 and 111 poles, which will not be the case for other crystals. But the spherical coordinate grid can not satisfy having points in the corners of the fundamental zone either, so in that respect the grids are equally valid.

I propose I throw my grid into its own separate function and I will utilize your cropping functionality - I'll work from a branch branched from this PR. Then I could make a PR and you can review it/we can discuss it further.

[din14970]

Ignore what I said about platonic polyhedra, the fairest sampling of a sphere can be achieved with a the vertices of a geodesic : https://en.m.wikipedia.org/wiki/Geodesic_polyhedron

[din14970]

I didn't realize the connection before, but essentially this is the problem of meshing the sphere, and there's multiple ways to do it each with their advantages drawbacks. I found the following article about 4 ways to do it: https://medium.com/game-dev-daily/four-ways-to-create-a-mesh-for-a-sphere-d7956b825db4. The currently implemented grid corresponds to the UV sphere, my methods correspond to the normalized cube (linear) and spherified cube (scaled). The icosahedral mesh is the most "fair" but also the most difficult to implement and offers the least control on "resolution". Maybe it would be nice to have all possible meshes implemented and offer the user the choice. The icosahedral mesh could maybe be implemented using this package: https://pypi.org/project/meshzoo/. I will try to do it in my PR. Once you have the mesh of the full sphere, cropping to the fundamental triangle is the same for all of them.
So maybe what I suggest is having different functions to generate each type of mesh/grid, and then have an extra argument e.g. mesh in get_beam_directions_grid with which the user can choose the type of grid he would like.

[pc494]

Moving this discussion to an enhancement issue as it will involve a function signature change. I think I would still like the code currently in this PR to be included for a patch fix, as it does still constitute a "bugfix" for the previous, broken implementation, even if we get rid of it in 0.4.0

edit: for clarity

[hakonanes]

Regarding the IPF corners, I found that the full pole figure, plotted according to the procedure mentioned in #128 (comment), wasn't covered in the cubic, orthorhombic and monoclinic case using MTEX' corners. I therefore kept the existing corners (written by @shogas) for these systems, which actually covered the full pole figure.

Don't know how satisfactory that is... Would very much like to see the literature @shogas used for the corners.

[pc494]

I'm going to delay merging this until I've had time to have another proper look. There's no rush.

[hakonanes]

Okay, feel free to change the corners again!

[pc494]

I think this, even as just syntax improvements, is worth having in