EBSD data with hex-grid issues

[SteffenBrinckmann]

Hey,
some EBSD scan in a hex-grid not a rectangular grind. Opening the file, editing work great. But as soon as plotting starts, the remapping fails.
Example file: https://drive.google.com/file/d/1BQEGyIZWoIUiRbsJJuJv1G399otYiRRM/view?usp=sharing

I can brute force map the initial data into a rectangular grid, but then all editing and saving will fail. It would be best to do that during plotting only

What is the suggested path?
Best, Steffen

[hakonanes]

Hi @SteffenBrinckmann!

Could you provide the code you run and the error message? I assume the error relates to the coordinate arrays.

[hakonanes]

By the way, CrystalMap is made for square grids. If you need CrystalMap to support hexagonal grids, I could review any PRs towards that goal.

[SteffenBrinckmann]

Here an example code for the failing case:

import matplotlib.pyplot as plt
from orix.io import load
from orix import plot

ebsd = load('EBSD.ang') #see download
fig = plt.figure()
ax = fig.add_subplot(projection="plot_map")
im = ax.plot_map(ebsd)  #wraps matplotlib.axes.Axes.imshow
plt.show()

For interpolation into a rectangular grid

import numpy as np
from scipy.interpolate import griddata
stepSize = ebsd.x[1]-ebsd.x[0]
xMax = np.max(ebsd.x)
xMin = np.min(ebsd.x)
yMax = np.max(ebsd.y)
yMin = np.min(ebsd.y)
widthPixel  = int((xMax-xMin)/stepSize)
ratio       = (xMax-xMin)/ (yMax-yMin)
heightPixel = int(widthPixel/ratio)
xAxis = np.linspace(xMin, xMax, widthPixel)
yAxis = np.linspace(yMin, yMax, heightPixel)
newX, newY = np.meshgrid(xAxis, yAxis)

oldPoints    = np.vstack(  (ebsd.x, ebsd.y) ).T
oldRotations = ebsd.rotations.to_euler()
newEuler0 = griddata( oldPoints, oldRotations[:,0], (newX,newY), 'linear').flatten()
newEuler1 = griddata( oldPoints, oldRotations[:,1], (newX,newY), 'linear').flatten()
newEuler2 = griddata( oldPoints, oldRotations[:,2], (newX,newY), 'linear').flatten()
newPhaseID= griddata( oldPoints, ebsd.phase_id,     (newX,newY), 'nearest').flatten()
newIQ     = griddata( oldPoints, ebsd.unknown1,     (newX,newY), 'linear').flatten()
newDP     = griddata( oldPoints, ebsd.unknown2,     (newX,newY), 'linear').flatten()
newX      = newX.flatten()
newY      = newY.flatten()

#Assemble new data set
from orix.crystal_map import CrystalMap
from orix.quaternion.rotation import Rotation
eulerAngles = np.column_stack((newEuler0, newEuler1, newEuler2))
rotations = Rotation.from_euler(eulerAngles)
properties = {"iq": newIQ, "dp": newDP}
ebsd2 = CrystalMap(rotations=rotations,
                 phase_id=newPhaseID,x=newX, y=newY,
                 phase_list=ebsd.phases, prop=properties)

[hakonanes]

Thanks, could you also print the error message you get?

[SteffenBrinckmann]

Error message:

Traceback (most recent call last):
  File "solution.py", line 15, in <module>
    im = ax.plot_map(ebsd)  #wraps matplotlib.axes.Axes.imshow
  File "/home/sbrinckm/FZJ/SourceCode/orix/orix/plot/crystal_map_plot.py", line 166, in plot_map
    phase_id = crystal_map.get_map_data("phase_id")
  File "/home/sbrinckm/FZJ/SourceCode/orix/orix/crystal_map/crystal_map.py", line 744, in get_map_data
    array[self.is_in_data] = data
IndexError: boolean index did not match indexed array along dimension 0; dimension is 47817 but corresponding boolean dimension is 23909

[hakonanes]

Hm, seems like some array in CrystalMap and the array to keep track of which points are "in the data", like only indexed points or points with a certain phase id, have different sizes. is_in_data has about half the number of elements as array (some property, like phase_id).

Will look into this. I'm not 100% sure it is a bug, it might be that some arrays are not set correct when reading the TSL orientation data from the .ang file with an hexagonal grid, because the .ang reader has only been tested on a square grid.

[hakonanes]

Just found this plotting of hexagonal grids with Matplotlib from pyEBSD (https://github.com/arthursn/pyebsd/blob/master/pyebsd/ebsd/plotting.py#L1157), might be of interest to look, in addition to your conversion, at if extending CrystalMap for hexagonal grids.

[SteffenBrinckmann]

I'm not sure, I would do it that complicated / slow way.
As far as I see:

• they plot filled hexagons and calculate the verticies for each. So if you have 1e6 datapoints, you have to draw 6e6 verticies, many of which are visually not able to distinguish and slow. That is what pyebsd is calling hex-tiling.
• I would interpolate into a regular grid and plot that as it is faster in plotting. That is called rect-tiling
• In a rectangular grid, you also do not draw rectangles but points.
  But that is only my 2ct.

[hakonanes]

they plot filled hexagons and calculate the verticies for each. So if you have 1e6 datapoints, you have to draw 6e6
verticies, many of which are visually not able to distinguish and slow. That is what pyebsd is calling hex-tiling.

I agree, that doesn't sound feasible!

I would interpolate into a regular grid and plot that as it is faster in plotting. That is called rect-tiling

Yes, I'm starting to think that we should have some CrystalMap class property for keeping track of grid type (rectangular or hexagonal), and at least support this interpolation in CrystalMap.get_map_data(). That method needs updating anyway. We should return a NotImplemented warning (or error) whenever a CrystalMap operation doesn't support an hexagonal grid, and then whoever has the time can start to implement the operations... I haven't thought much about this, but I hope we can start from this without too much effort.

[pc494]

Sorry to jump in at this late stage. Am I correct in thinking that a hex2rect algorithm would need to interpolate on all properties of the CrystalMap? I think it might be a challenge to write a one size fits all solution for this.

Examples:
How do you interpolate a discrete value (such as phase)?
How do you interpolate Orientations that are far apart from one another (and may have symmetry relations)?

[hakonanes]

Sorry to jump in at this late stage. Am I correct in thinking that a hex2rect algorithm would need to interpolate on all properties of the CrystalMap? I think it might be a challenge to write a one size fits all solution for this.

Yes, I would think so, and I totally agree that this should be carefully handled.

To address hexagonal grids, we could start by raising a warning when loading from file (TSL stores grid type in their .ang file), or raise a NotImplementedError. The .ang reader is not tested for hexagonal grids, so I have no idea how the resulting CrystalMap looks. But it is apparently possible to read hexagonal .ang files into a CrystalMap, given this issue's topic.

I would rather it be possible to read the data, but explicitly warn that this isn't supported yet, and then slowly start to support it. Do you agree @pc494?

[pc494]

This sounds like a good compromise while we consider the most suitable way to handle such data (slash await time to work on such a project)

[hakonanes]

@IMBalENce brought square to hexagonal grid conversion in MTEX, https://mtex-toolbox.github.io/EBSDGrid.html, to my attention. This implementation might be useful to look at in addition to the one mentioned above #214 (comment).

[hakonanes]

Short term solution here is to warn when orientations in an .ang file are stored in an hexagonal grid. The warning should state that the CrystalMap class isn't made for hexagonal grids, and thus the behaviour of its properties and methods is undefined. I think this should be added to v0.6.1, and then those people that work with hexagonal grids can contribute to make orix support it.

[pc494]

I'm happy to do this

[hakonanes]

Was contacted via e-mail by a potential user who has data acquired in a hexagonal grid. Just mentioning it here to make issue readers aware of interest in this feature.

[argerlt]

Just a comment on this:

they plot filled hexagons and calculate the verticies for each. So if you have 1e6 datapoints, you have to draw 6e6 verticies, many of which are visually not able to distinguish and slow.

This is actually what plt.imshow does in the background, but with square grids instead of hexes. It feels slow until you realize it's just easier ray tracing.

Additionally, This is how MTEX creates all EBSD plots, be they square, hex, or irregular. First it uses a program called SpatialDecomposition.m to convert pixel coordinates into lists of nodes, vertices, and faces, then plots the EBSD scans as 2D planar graphs.

Problem is, you can't quickly do these sorts of plots in Matplotlib. Instead, you need modules that use openGL or something equivalent, like pyvista, glumpy, or pyvis. Also, there is no pythonic equivalent to SpatialDecomposition, though scipy.spatial.Delaunay is close.

I will add, this is ONE way to support hex-grid ebsd, but certainly not the ONLY way. However, this methodology is also what I was talking about in #151 as a way to quickly make region maps, so creating a class for node/vertex/face maps would have multiple benefits.

Also, to give some context, consider this glumpy example that plots 1E6 spheres moving at 60 fps:
https://github.com/glumpy/glumpy/blob/master/examples/gloo-cloud.py