IPF color maps

[din14970]

Since the addition of stereographic plots, orix has become a lot more useful also for visualisation. A really nice addition to the visualisation capabilities would be quick IPF maps: I want to be able to iterate on my template matching parameters through feedback I get from a visual, and the most intuitive visualisation are IPF color maps. I have worked out some stuff for my specific case of a cubic crystal, but throwing this as-is into orix is probably not ideal. Nevertheless, I thought I'd throw what I have in here to kickstart a discussion.

In principle mapping vectors to a color is not difficult, but I have yet to figure out how those stereographic triangle colormaps are actually implemented. For the cubic case I came up with something which looks reasonable and close to the map in mtex:

def get_ipf_color(vectors):
    """Map vectors to a color"""
    # the columns of this matrix should map to red, green blue respectively
    color_corners = np.array([[0, 1/np.sqrt(2), 1/np.sqrt(3)],
                              [0, 0, 1/np.sqrt(3)],
                              [1, 1/np.sqrt(2), 1/np.sqrt(3)]])
    color_mapper = np.linalg.inv(color_corners)

    # a bit of wrangling
    data_sol = vectors.data
    flattened = data_sol.reshape(np.product(data_sol.shape[:-1]), 3).T
    rgb_mapped = np.dot(color_mapper, flattened)
    rgb_mapped = np.abs(rgb_mapped / rgb_mapped.max(axis=0)).T
    rgb_mapped = rgb_mapped.reshape(data_sol.shape)**0.5
    return rgb_mapped

On a cubic grid:

The MTEX color map:

It's not a perfect match, but it does the job.

I can then use this function to map my orientation mapping results to a color image:

from orix.quaternion.rotation import Rotation
from orix.vector.vector3d import Vector3d

# The matrix
solution_vectors_1 = Rotation.from_euler(np.deg2rad(result_steel["austenite"]["orientation"][:,:,0,:]))*Vector3d.zvector()
# The particle
solution_vectors_2 = Rotation.from_euler(np.deg2rad(result_gph["g"]["orientation"][:,:,0,:]))*Vector3d.zvector()

plt.imshow(get_ipf_color(solution_vectors_1))
plt.imshow(get_ipf_color(solution_vectors_2))

I think it would be cool to get something like this working out of the box on orix. Any ideas on how to go about it?

[hakonanes]

This is a very good starting point for inverse pole figure (IPF) coloring, @din14970, thanks!

We have descriptions of all symmetry operations in the 32 crystallographic point groups in all seven crystal systems in the orix.quaternion.symmetry module. We should therefore work towards pole figure and IPF functionality, including IPF coloring, supporting all these groups. So I agree that we shouldn't throw your code in as-is. In my view, with orix, we should implement tools, not routines. The latter are implemented in packages like pyxem and kikuchipy.

With that said, I'm almost done with the Miller class in https://github.com/pyxem/orix/pull/164 for all point groups. The natural continuation is to implement pole figure plotting of orientations, and also restricting (zooming) the stereographic projection to within three vectors (perhaps even just two if they their polar angle is pi / 2), which we then can use to get the fundamental sector/triangle for every point group. This isn't far off.

@din14970, could you paste here the full code you use to plot your IPF coloring? And since we in orix plot x to east, shouldn't your vectors in get_ipf_color() be [001], [101], [111]?

You should be able to import the classes Rotation and Vector3d like so:

from orix.quaternion import Rotation
from orix.vector import Vector3d

Side-note: Since implementing Vector3d.scatter() and Vector3d.draw_circle() turned out so simple, I have quite immediate plans to implement CrystalMap.plot(), to make this plotting MUCH easier.

[din14970]

And since we in orix plot x to east, shouldn't your vectors in get_ipf_color() be [001], [101], [111]?

Basically in the first step I am relying on a linear transformation of:

[001] (first vector) --> [001] (normalized) --> [100] (red)
[101] (second vector) --> [1/sqrt(2) 0 1/sqrt(2)] --> [010] (green)
[111] (third vector  --> [1/sqrt(3) 1/sqrt(3) 1/sqrt(3)] --> [001] (blue)

I figured normalizing the vectors would be the most fair, if one does not normalize, and use

color_corners = np.array([[0, 1, 1],
                                           [0, 0, 1],
                                           [1, 1, 1]])

Then the map looks like

So basically by adjusting the length of the vectors one can stretch the color wheel to different colors.

In the second bit under a bit of wrangling I am ensuring that in all the colors, at least one of RGB has the value of 1 and the others less. The square root serves to spread out the distribution a bit more, otherwise most orientations are primarily red, green and blue. With the root removed:

You should be able to import the classes Rotation and Vector3d like so

The reason I often do it like this is because my tab-autocomplete in a jupyter notebook finds the things in the actual file or folder, but not those that are imported via __init__.py. So it is sometimes more convenient to do the less convenient thing ;).

Full code for IPF coloring

For the stereogram:

import numpy as np
from orix.quaternion import Rotation
from orix.vector import Vector3d
from orix.projections import StereographicProjection
from diffsims.generators.rotation_list_generators import get_beam_directions_grid

# grid of euler angles
resolution = 1 
grid_cub = get_beam_directions_grid("cubic", resolution, mesh="spherified_cube_edge")

# grid of vectors
grid_stereo = Rotation.from_euler(np.deg2rad(grid_cub))*Vector3d.zvector()

sp = StereographicProjection()

fig, ax = plt.subplots(figsize=(8,8))
ax.set_aspect("equal")
ax.scatter(*sp.vector2xy(grid_stereo), s = 50, c=get_ipf_color(grid_stereo.unit))  # use get_ipf_color for color of points

Sidenote: should vector2xy not be a static method?

For the IPF maps it's a similar procedure

result = index_dataset_with_template_rotation(data, diff_lib)
rotations = result["austenite"]["orientation"][:,:,0,:]  # euler angles as (x, y, 3) array. The 0 is to select the first solution.
vecs_z = Rotation.from_euler(np.deg2rad(rotations))*Vector3d.zvector()  # z vectors in a 2D grid to be plot
plt.imshow(get_ipf_color(vecs_z))  # convert vectors to RGB image

[hakonanes]

Being able to easily stretch the color key is very important, and should be a feature!

~What I was saying with the crystal directions [uvw] was that your colored triangle has corners [001], [101], and [111], at least in our stereographic plot, right? If you use [011], you would get another triangle, basically the next if you rotate counter-clockwise.~ I totally misunderstood, mistaking RGB values for crystal directions...

The gridding functions in diffsims are on S2, right? I think these should be in the orix.sampling module, what do you think?

Sidenote: should vector2xy not be a static method?

No, we keep track of whether vectors are in the current spherical region (for now upper or lower hemisphere) with the StereographicProjection.region property: https://github.com/pyxem/orix/blob/master/orix/projections/stereographic.py#L49. This way we can in the future specify which pole, i.e. orix.vector.SphericalRegion, of the unit sphere we look down, not just northern or southern hemispheres.

---

When we implement a IPFKey class, I think it should have some important properties/features like:

• direction: with respect to which sample reference direction should the vectors be colored? Basically your vector multiplication.
• color key type, say TSL or something like MTEX has (http://dx.doi.org/10.1107/S1600576716012942)
• should be initialized with a point group (different triangles)
• should take in orientations (or rotations), and return an RGB array

@din14970, what do you think of this?

[din14970]

The gridding functions in diffsims are on S2, right? I think these should be in the orix.sampling module, what do you think?

Yes, I have the constraint that phi1 = 0 in my grid function and so am basically sampling a 3D sphere surface. Sometimes it's hard to see the scope of the different packages, but yes I might agree with you that sampling orientations might fit better in orix. One issue is that I do rely on the orientations being in euler angle format as a numpy array, not an orientation object, for passing it into the template matching algorithm because template matching attempts to find phi1. If I turn it into an orientation and back into euler angles, the first euler angle is often no longer 0 because there are often multiple ways to arrive at the same orientation. How these diffraction libraries are generated with diffsims is quite messy and convoluted in my opinion, this might be a point to work on.

@din14970, what do you think of this?

Yep absolutely right, input is orientations/rotations and output is RGB. One should indeed have direction so that IPF-X and Y are easy, or any arbitrary vector really. These things are simple. The challenge is indeed how to map the colors and how to do it for all point groups, and since there appears to be a paper about it it's not so trivial. For no symmetry it's trivial, then it's basically mapping the HSV color cylinder to the sphere. But with the triangles I'm confused about how to properly do it.

[pc494]

The gridding functions in diffsims are on S2, right? I think these should be in the orix.sampling module, what do you think?

Yes, I have the constraint that phi1 = 0 in my grid function and so am basically sampling a 3D sphere surface. Sometimes it's hard to see the scope of the different packages, but yes I might agree with you that sampling orientations might fit better in orix. One issue is that I do rely on the orientations being in euler angle format as a numpy array, not an orientation object, for passing it into the template matching algorithm because template matching attempts to find phi1. If I turn it into an orientation and back into euler angles, the first euler angle is often no longer 0 because there are often multiple ways to arrive at the same orientation. How these diffraction libraries are generated with diffsims is quite messy and convoluted in my opinion, this might be a point to work on.

I would prefer that functionality to stay in diffsims for now. orix.sampling currently contains functionalities that do "proper samplings". The diffsims "sampling" is based on exploiting a specific property of SED patterns (that pattern simulations are related by an image rotation if phi1 = 0) and so keeping that stuff sectioned off made sense. However, if our colouring scheme relies on the same property (that [001] is red, regardless of how much we've spun it around) then we should probably shift it over. Worth being careful about what we name things as we do that though.

[hakonanes]

orix.sampling currently contains functionalities that do "proper samplings". The diffsims "sampling" is based on exploiting a specific property of SED patterns (that pattern simulations are related by an image rotation if phi1 = 0) and so keeping that stuff sectioned off made sense.

What I mean is that the spherical gridding functions should be upstreamed to orix, since this is very nice when e.g. generating vectors to create color mappings as @din14970 has done above or like I've done here for point group m-3m (note that Matplotlib in figure.savefig() makes the scatter points quite large and overlapping, so the color triangle looks semi-continuous)

I created the delineating triangle with the new Vector3d.get_circle() and StereographicPlot.draw_circle() and restricting the azimuthal and polar angles using the orix.vector.SphericalRegion defined by the [1, 1, 0], [1, 1, 1] and [0, 0, 1] crystal plane normals.

I used your color mapping function for point group m-3m, @din14970, and got these results:

The pole figures show basically, in pseudo code, rotations * vector.Miller(hkl=[1, 0, 1], phase=m-3m) etc. What is needed next is a projection to the fundamental sector, which is the SphericalRegion delineated in the top figure.

[din14970]

Looks cool! Have you compared the colors to what you get with MTEX?

[hakonanes]

No, I haven't compared to colors from MTEX yet, because we must take great care that symmetrically equivalent rotations and directions are handled correctly... This isn't far off, but not there yet.

My goal is to get #164 done first, then move on to plotting symmetry elements from point groups, basically a Symmetry.plot(), where we also can delineate the fundamental sector, as shown in the large stereographic plot above. Will open an issue/PR about this soon. I need this to get a visual on various symmetries and fundamental sectors, before creating an IPFKey class to color orientations.

[hakonanes]

Having looked at this some more, both what is done in MTEX and in the above referenced paper by Nolze and Hielscher, we can get these colours for point group m-3m for an Ni dataset with orix comparing to MTEX' TSL color key

There are subtle differences which I cannot explain, but that is something we can iron out as we start to get PRs in with this functionality.

Computing ~orientations * x/y/zvector, projecting them to the fundamental region, and then coloring the returned vectors, we can get this figure

There are a lot of tools (classes/functions/methods) required to obtain this result, mainly:

• The Miller class in #164
• A description of each point group's fundamental sector on the unit sphere (I didn't get the expected vertices from Symmetry.fundamental_sector() already implemented, e.g. point group m-3m returned an empty vector set), with vertices and a center. In the stereographic projections above, I plotted the great circles of the vectors describing the fundamental sector of point group m-3m. By taking the vector cross products we get the vertices, which are also plotted.
• Projection of vectors to the fundamental sector of a point group.
• Coloring of the projected vectors:
  • Calculation of the projected vectors' azimuth and polar angles relative to the fundamental sector center (refilling the stereographic projection with the -center as stereographic projection pole)
  • Fitting the HSV (hue, saturation, value) color map to each vector relative to the sector center

What's nice is that orix' vector and symmetry classes have almost all we need, apart from the Symmetry.fundamental_sector(), so it's not that many new lines of code.

[hakonanes]

The blue dots in the above projection are the vectors before projecting to the fundamental sector. The figure basically shows the inverse pole figure.

[din14970]

Looks cool!

There are subtle differences which I cannot explain, but that is something we can iron out as we start to get PRs in with this functionality.

That's because my way of mapping the colors is done based on the corners and linear algebra, whereas it seems the MTEX way is based on the centroid of the triangle which is mapped to white, and then the HSV color wheel is somehow stretched inside the triangle. I don't have a clue how to implement it the MTEX way but consistency would be nice.

Coloring of the projected vectors: Calculation of the projected vectors' azimuth and polar angles relative to the fundamental sector center (refilling the stereographic projection with the -center as stereographic projection pole). Fitting the HSV (hue, saturation, value) color map to each vector relative to the sector center

I think this will be much more tricky than expected and I'm not sure whether polar coordinates will be an appropriate parametrization. The (101) and (111) poles are not on the same latitude, so in m-3m ensuring that (101) is pure green and (111) is pure blue with a pure hue gradient in between I think will be very tricky when using a polar coordinate map. Honestly what would you think of just posting on the MTEX forums to ask how they handle the colors?

[hakonanes]

That's because my way of mapping the colors is done based on the corners and linear algebra, whereas it seems the MTEX way is based on the centroid of the triangle which is mapped to white, and then the HSV color wheel is somehow stretched inside the triangle. I don't have a clue how to implement it the MTEX way but consistency would be nice.

sorry, i should have said this, but i didn't use your mapping @din14970, but what is explained in the referenced paper and how it is done in mtex. that also goes for designating hsv colours to the polar coordinates.

what im saying is that i have the functionality as simple functions, but i want them structured into classes and more general functions so we can build upon this!

[din14970]

sorry, i should have said this, but i didn't use your mapping @din14970, but what is explained in the referenced paper and how it is done in mtex. that also goes for designating hsv colours to the polar coordinates.

Ah sorry about that, so you already figured it out then, awesome!

[hakonanes]

I should mention that after projecting all rotated vectors to the fundamental region, your coloring isn't too far off MTEX' TSL coloring, @din14970.

[dnjohnstone]

Just to say - this is awesome!

[ravipurohit1991]

Hello, can we plot the IPF color codes with the recent version of Orix (however I did not find the relevant class/functions in the master) or this is still under development ? I would very much appreciate this feature in the Orix library. Thank you.

[hakonanes]

Hi @ravipurohit1991 ! Unfortunately, we can't do that yet, but we're working on it. As I mention in this issue, there are a few things needed before we can do this: https://github.com/pyxem/orix/issues/166#issuecomment-808155612

Thanks for showing interest!

[harripj]

This would be a great feature! Just to add that it should be possible to use plt.pcolorfast to generate a coloured mesh plot rather than coloured points if the orientations can be arranged into a square grid. Here is something I threw together:

# have some orientations (rot, scipy in this case) from a grid (of size side x side) that represent the standard triangle
xyz = rot.apply((0, 0, 1)).reshape(side, side, 3)
# split into components
x, y, z = xyz[..., 0], xyz[..., 1], xyz[..., 2]

# get grid center points for colour evaluation
x_mid = np.diff(np.diff(x, axis=0), axis=1) + x[:-1, :-1]
y_mid = np.diff(np.diff(y, axis=0), axis=1) + y[:-1, :-1]
z_mid = np.diff(np.diff(z, axis=0), axis=1) + z[:-1, :-1]
xyz_mid = np.column_stack((x_mid.ravel(), y_mid.ravel(), z_mid.ravel()))

c = get_ipf_color(xyz_mid).reshape(*x_mid.shape, 3)

fig, ax = plt.subplots()
# project and show colours
ax.pcolorfast(x / (1 + z), y / (1 + z), c.clip(0, 1))

NB. I have clearly introduced some skew in the orientation grid when trying to create it. If someone knows of a better way that would be great!

[hakonanes]

that's a great tip, @harripj! the mesh method should work fine with the StereographicAxes. ill have a look at it.

[din14970]

alternatively there is also scipy.interpolate.griddata to plot a smooth image instead of individual points

from scipy.interpolate import griddata
from diffsims.generators.rotation_list_generators import get_beam_directions_grid
import matplotlib.colors as mcolors
from orix.projections import StereographicProjection

resolution = 1
grid_cub = get_beam_directions_grid("cubic", resolution, mesh="spherified_cube_edge")

# project a set of vectors to the fundamental zone of the cubic crystal
def to_fundamental(data_sol):
    data_sol = np.abs(data_sol)
    data_sol = np.sort(data_sol, axis=-1)
    column = data_sol[...,0].copy()
    data_sol[..., 0] = data_sol[...,1]
    data_sol[..., 1] = column
    return data_sol

# get the color corresponding to a set of vectors
def get_ipf_color(vectors):
    # the following column vectors should map onto R [100], G [010], B[001], i.e. the identity. So the inverse of 
    # this matrix maps the beam directions onto the right color vector
    color_corners = np.array([[0, 1, 1],
                              [0, 0, 1],
                              [1, 1, 1]])
    color_mapper = np.linalg.inv(color_corners)

    # a bit of wrangling
    data_sol = to_fundamental(vectors.data)
    flattened = data_sol.reshape(np.product(data_sol.shape[:-1]), 3).T
    rgb_mapped = np.dot(color_mapper, flattened)
    rgb_mapped = np.abs(rgb_mapped / rgb_mapped.max(axis=0)).T
    rgb_mapped = rgb_mapped.reshape(data_sol.shape)
    return rgb_mapped

# a helper function for turning the euler angle grid to a grid of points in the stereographic projection
def grid_to_xy(grid):
    from orix.quaternion.rotation import Rotation
    from orix.vector.vector3d import Vector3d
    from orix.projections import StereographicProjection
    s = StereographicProjection(pole=-1)
    rotations_regular =  Rotation.from_euler(np.deg2rad(grid))
    rot_reg_test = rotations_regular*Vector3d.zvector()
    x, y = s.vector2xy(rot_reg_test)
    return x, y

# getting the z vector corresponding to an orientation
def ori_to_vec(eulers):
    from orix.quaternion.rotation import Rotation
    from orix.vector.vector3d import Vector3d
    rotations_regular =  Rotation.from_euler(np.deg2rad(eulers))
    return rotations_regular*Vector3d.zvector()

# xy coordinates of the points in the stereographic projection
xy = np.array(grid_to_xy(grid_cub)).T
colors = get_ipf_color(ori_to_vec(grid_cub))
reds = colors[:, 0]
greens = colors[:, 1]
blues = colors[:, 2]

# here we interpolate between these points on a square grid. We do it for red, green and blue separately 
sampling=0.001
gridx, gridy = np.mgrid[-0.05:0.42:sampling, -0.05:0.45:sampling]
t_rd = griddata(xy, reds, (gridy, gridx), method="linear")
t_gn = griddata(xy, greens, (gridy, gridx), method="linear")
t_bl = griddata(xy, blues, (gridy, gridx), method="linear")
t_alpha = np.invert(np.isnan(t_rd))  # everything outside the triangle is black and transparant
t_rd[np.isnan(t_rd)] = 0
t_bl[np.isnan(t_bl)] = 0
t_gn[np.isnan(t_gn)] = 0
triangle = np.stack([t_rd, t_gn, t_bl, t_alpha], axis=-1)

fig, ax1 = plt.subplots()
ax1.imshow(triangle)
ax1.invert_yaxis()

[harripj]

That looks great!

[hakonanes]

Indeed, that looks very nice. We should try to use these functions when plotting the IPF color map reference.

[hakonanes]

See the follow-up PR is #235.