How to plot a phase map

[EffettiCollaterali]

Dear all
I'm working with Kikuchipy with the aim of distinguish Fe4N from Ferrite, and I would like to use, in a second moment, the dictionary indexing refinement to do that. At the moment I'm stuck on how to plot a phase map after the indexing to see at least if two phases are somewhat distinguished. Is there also a method to print the percentage of indexed points for each phase?

Thanks in advance,

Andrea

The following is a test piece of code that I'm using (the reflectors part was taken from this forum).

phase_list = PhaseList(
    names=["Fe", "Al"],
    space_groups=[215, 229],
    structures=[
        Structure(
            atoms=[Atom("fe", [0] * 3)],
            lattice=Lattice(0.379, 0.379, 0.379, 90, 90, 90)
        ),
        Structure(
            atoms=[Atom("al", [0] * 3)],
            lattice=Lattice(0.286, 0.286, 0.286, 90, 90, 90)
        ),
    ]
)

reflectors = []
for indx, ph in phase_list:
    if indx != -1:
        lattice = ph.structure.lattice
        lattice.setLatPar(lattice.a * 10, lattice.b * 10, lattice.c * 10)

        rlv = ReciprocalLatticeVector.from_min_dspacing(ph, 0.7)
        rlv.sanitise_phase()
        rlv = rlv.unique(use_symmetry=True)
        rlv.calculate_structure_factor("lobato")
        structure_factor = abs(rlv.structure_factor)
        order = np.argsort(structure_factor)[::-1]
        rlv = rlv[order[0:7]]    # I workaround: slice the array because it contains values giving ValueError: negative dimensions are not allowed
        rlv = rlv.symmetrise()
        rlv.print_table()
        reflectors.append(rlv)

indexer = det.get_indexer(
    phase_list=phase_list, 
    reflectors=reflectors, 
    nBands = 10,
    tSigma = 2,
    rSigma = 2
)

print(indexer.vendor)
print(indexer.sampleTilt)
print(indexer.camElev)
print(indexer.PC)

print(indexer.phaselist[0].latticeparameter)
print(indexer.phaselist[0].polefamilies)
print(indexer.phaselist[1].latticeparameter)
print(indexer.phaselist[1].polefamilies)

det = s_grid.hough_indexing_optimize_pc(
    pc0=[0.42, 0.22, 0.50],  # Initial guess based on previous experiments
    indexer=indexer,
    batch=True,
    method="PSO",
    search_limit=0.05,
)

# Print mean and standard deviation
print(det.pc_flattened.mean(axis=0))
print(det.pc_flattened.std(0))

det.plot_pc("scatter", s=50, annotate=True)
det.pc = det.pc_average
det.plot(pattern=s_grid.inav[0, 0].data)

indexer = det.get_indexer(
    phase_list=phase_list, 
    reflectors=reflectors, 
    nBands = 10,
    tSigma = 2,
    rSigma = 2
)
indexer.PC

xmap = s.hough_indexing(phase_list=phase_list, indexer=indexer, verbose=2)

[hakonanes]

Hi Andrea!

distinguish Fe4N from Ferrite

Given these phase descriptions:

• Fe4N (source): space group 221 (Pm-3m), a = 0.376 nm, N in the corners and Fe on the edge centers and in the center
• Ferrite: space group 229 (Im-3m), a = 0.287 nm, BCC

I don't think we can distinguish between these two phases using Hough indexing. This is because they have the same Laue group, m-3m, and thus the same interplanar angles between high-symmetry planes. See e.g. this application note from EDAX. There, they had to use energy dispersive X-ray spectroscopy to distinguish the nitride from ferrite.

You may, on the other hand, be able to differentiate between the two based on dynamical simulations! I haven't seen simulations for Fe4N myself, but given the differing crystal structures (lattice parameters, the nitrogen atoms), it may be possible.

---

Regarding your technical questions, EBSD.hough_indexing() returns a CrystalMap. You can get the phase fractions by just inspecting the xmap instance. For a crystal map xmap from super-duplex stainless steel (iq = image quality, dp = dot product values from dictionary indexing with EMsoft):

>>> xmap
Phase   Orientations       Name  Space group  Point group  Proper point group       Color
    1   5657 (48.4%)  austenite         None          432                 432    tab:blue
    2   6043 (51.6%)    ferrite         None          432                 432  tab:orange
Properties: iq, dp
Scan unit: um

You can then plot the map using CrystalMap.plot()

>>> xmap.plot()

I hope this helps!

[EffettiCollaterali]

Dear Håkon,
thank you very much for the help. I was able to plot the data. The grains of my material are small (200 nm), at the moment I'm able to distinguish austenite and martensite grains only using transmission EBSD. Is that a problem for the indexing method to use the pattern obtained in transmission? The sample is tilted negatively 20 degrees in a Bruker system. In the same way, I would like to rely on transmission data also for the Fe4N grain discrimination.

Best,
Andrea
(FBK, Trento, Italy)

[hakonanes]

Good question. There shouldn't be, but I haven't tried it. If you could share a single pattern from austenite with estimated projection parameters, I could have a go. Is that possible?

[EffettiCollaterali]

Dear Håkon,
Do you mean the PCs from Kikuchipy or the Bruker software?
This is one of the EBSD files in transmission.
https://drive.google.com/file/d/19TgfKldvjJt10nsdgqoe62smdm8B88bX/view?usp=drive_link
This is the phase map, that I obtained with the Bruker software, where Fe4N is rendered in blue. I was surprised that a whole grain was indexed as Fe4N. The indexing of that grain as Fe4N is consistent even increasing the number of bands up to 10.

I tried to use Kikuchipy for the same data, however, I am still far from getting similar results.

The following is the code that I'm currently using.

import matplotlib.pyplot as plt
import numpy as np
import pyvista
import skimage.exposure as ske
import skimage.transform as skt
import hyperspy.api as hs
import kikuchipy as kp
from orix import io, plot, quaternion, vector
from orix.quaternion import Orientation, Rotation
from orix.vector import Vector3d
from diffpy.structure import Atom, Lattice, Structure
from diffsims.crystallography import ReciprocalLatticeVector
import kikuchipy as kp
from orix import plot
from orix.crystal_map import Phase, PhaseList
from orix.vector import Vector3d
import hyperspy.api as hs

#Size of the images

plt.rcParams["figure.figsize"] = (15, 5)

plt.rcParams.update({"font.size": 15, "lines.markersize": 15})

#Load the data

s=kp.load("Steel.hdf5")
#s = kp.data.nickel_ebsd_large(allow_download=True)

print(s.static_background)


#Plot the data

s.plot()
plt.show()


#Downsize the patterns

s2 = s.downsample(2, inplace=False)

s = s2


#Remove the dynamic background

s.remove_dynamic_background()
s.plot()
plt.show()


#Remove the static background

static_bg = s.mean(axis=(0, 1))  # (0, 1) are the navigation axes
static_bg.rescale_intensity(dtype_out=np.uint16)
s.remove_static_background(static_bg=static_bg.data)
s.plot()
plt.show()

# Average in a (3, 3) window with a Gaussian kernel with a standard
# deviation of 1
s.average_neighbour_patterns(window="gaussian", std=2)
s.plot()
plt.show()

iq1 = s.get_image_quality()


print(iq1.mean())

plt.imshow(iq1, cmap="gray")
plt.tight_layout()

dtype_orig = s.data.dtype
s.change_dtype("float64")

#n_components = 100
#s.decomposition(
#    algorithm="SVD",
#    output_dimension=n_components,
#    centre="signal",
#)

#s.change_dtype(dtype_orig)

#s.plot_decomposition_results()

#factors = s.learning_results.factors  # (n detector pixels, m components)
#sig_shape = s.axes_manager.signal_shape[::-1]
#loadings = s.learning_results.loadings  # (n patterns, m components)
#nav_shape = s.axes_manager.navigation_shape[::-1]

#fig, ax = plt.subplots(ncols=2, figsize=(10, 5))
#ax[0].imshow(loadings[:, 0].reshape(nav_shape), cmap="gray")
#ax[0].axis("off")
#ax[1].imshow(factors[:, 0].reshape(sig_shape), cmap="gray")
#ax[1].axis("off")
#fig.tight_layout(w_pad=0)

#_ = s.plot_explained_variance_ratio(n=n_components)

#fig, ax = plt.subplots(ncols=4, figsize=(15, 5))
#for i in range(4):
#    factor_idx = i + 90
#    factor = factors[:, factor_idx].reshape(sig_shape)
#    factor_iq = kp.pattern.get_image_quality(factor)
#    ax[i].imshow(factor, cmap="gray")
#    ax[i].set_title(f"#{factor_idx}, IQ = {np.around(factor_iq, 2)}")
#    ax[i].axis("off")
#fig.tight_layout()

#s2 = s.get_decomposition_model(components=90)

#iq2 = s2.get_image_quality()
#iq2.mean()


#fig, ax = plt.subplots(ncols=2, figsize=(15, 4))
#im0 = ax[0].imshow(iq1, cmap="gray")
#ax[0].axis("off")
#fig.colorbar(im0, ax=ax[0], pad=0.01, label="IQ before denoising")
#im1 = ax[1].imshow(iq2, cmap="gray")
#ax[1].axis("off")
#fig.colorbar(im1, ax=ax[1], pad=0.01, label="IQ after denoising")
#fig.tight_layout(w_pad=-10)

#hs.plot.plot_signals([s, s2], navigator=hs.signals.Signal2D(iq2))

#s = s2

#Set the detector

det = kp.detectors.EBSDDetector(
    shape=s.axes_manager.signal_shape[::-1],
    pc=[0.5, 0.4, 0.5],
    px_size=0.003,  # Microns
    #binning=4,
    tilt=7.7,
    sample_tilt=-18,
)

det.pc_tsl()

det.plot(pattern=s.inav[0, 0].data)


#Plot the Virtual BSE 

vbse_imager = kp.imaging.VirtualBSEImager(s)
print(vbse_imager.grid_shape)

maps_vbse_rgb = vbse_imager.get_rgb_image(r=(2, 1), g=(2, 2), b=(2, 3))

maps_vbse_rgb.plot()
plt.show()





#Plot the image quality

maps_iq = s.get_image_quality()
#s.xmap.plot(
#    maps_iq.ravel(),  # Must be 1D
#    cmap="gray",
#    colorbar=True,
#    colorbar_label="Image quality $Q$",
#    remove_padding=True,
#)
#plt.show()

# Extract data, also getting the grid positions
grid_shape = (6, 4)
s_grid, idx = s.extract_grid(grid_shape, return_indices=True)
s_grid

nav_shape = s.axes_manager.navigation_shape[::-1]

kp.draw.plot_pattern_positions_in_map(
    rc=idx.reshape(2, -1).T,  # Shape (n patterns, 2)
    roi_shape=nav_shape,  # Or maps_iq.shape
    roi_image=maps_iq,
)


phase_list = PhaseList(
    names=["Ferrite", "Austenite", "Fe4N"],
    space_groups=[229, 225, 221],
    structures=[
        Structure(
            atoms=[Atom("fe", [0] * 3)],
            lattice=Lattice(0.287, 0.287, 0.287, 90, 90, 90)
        ),
        Structure(
            atoms=[Atom("fe", [0] * 3)],
            lattice=Lattice(0.356, 0.356, 0.356, 90, 90, 90)
        ),
        Structure(
            atoms=[Atom("fe", [0] * 3), Atom("n", [0.5, 0.5, 0.5]), Atom("fe", [0, 0.5, 0.5])],
            lattice=Lattice(0.376, 0.376, 0.376, 90, 90, 90)
        ),
    ]
)

reflectors = []
for indx, ph in phase_list:
    if indx != -1:
        lattice = ph.structure.lattice
        lattice.setLatPar(lattice.a * 10, lattice.b * 10, lattice.c * 10)

        rlv = ReciprocalLatticeVector.from_min_dspacing(ph, 0.5)
        rlv.sanitise_phase()
        rlv = rlv.unique(use_symmetry=True)
        rlv.calculate_structure_factor("lobato")
        structure_factor = abs(rlv.structure_factor)
        order = np.argsort(structure_factor)[::-1]
        rlv = rlv[order[0:7]]    # I workaround: slice the array because it contains values giving ValueError: negative dimensions are not allowed
        rlv = rlv.symmetrise()
        rlv.print_table()
        reflectors.append(rlv)

indexer = det.get_indexer(
    phase_list=phase_list, 
    reflectors=reflectors, 
    nBands = 8,
    tSigma = 2.0,
    rSigma = 2.0
)

print(indexer.vendor)
print(indexer.sampleTilt)
print(indexer.camElev)
print(indexer.PC)

print(indexer.phaselist[0].latticeparameter)
print(indexer.phaselist[0].polefamilies)
print(indexer.phaselist[1].latticeparameter)
print(indexer.phaselist[1].polefamilies)

det = s_grid.hough_indexing_optimize_pc(
    pc0=[0.4591, 0.3679, 0.4611],  # Initial guess based on previous experiments
    indexer=indexer,
    batch=True,
    method="PSO",
    search_limit=0.05,
)

# Print mean and standard deviation
print(det.pc_flattened.mean(axis=0))
print(det.pc_flattened.std(0))

det.plot_pc("scatter", s=50, annotate=True)
det.pc = det.pc_average
det.plot(pattern=s_grid.inav[0, 0].data)

indexer = det.get_indexer(
    phase_list=phase_list, 
    reflectors=reflectors, 
    nT = 180,
    nBands = 8,
    tSigma = 2.0,
    rSigma = 2.0
)
indexer.PC

xmap = s.hough_indexing(phase_list=phase_list, indexer=indexer, verbose=2)

xmap

aspect_ratio = xmap.shape[1] / xmap.shape[0]
figsize = (8 * aspect_ratio, 4.5 * aspect_ratio)


fig, axes = plt.subplots(nrows=2, ncols=2, figsize=figsize, layout="tight")
for ax, to_plot in zip(axes.ravel(), ["pq", "cm", "fit", "nmatch"]):
    im = ax.imshow(xmap.get_map_data(to_plot))
    fig.colorbar(im, ax=ax, label=to_plot)
    ax.axis("off")


fig, axes = plt.subplots(ncols=2, nrows=2, figsize=(10, 5), layout="tight")
for ax, to_plot in zip(axes.ravel(), ["pq", "cm", "fit", "nmatch"]):
    ax.hist(xmap.prop[to_plot], bins=100)
    ax.set(xlabel=to_plot, ylabel="Frequency");


v_ipf = Vector3d.xvector()
sym = xmap.phases[0].point_group

ckey = plot.IPFColorKeyTSL(sym, v_ipf)
ckey

rgb_x = ckey.orientation2color(xmap.rotations)
fig = xmap.plot(rgb_x, overlay="cm", remove_padding=True, return_figure=True)

# Place color key in bottom right corner, coordinates are [left, bottom, width, height]
ax_ckey = fig.add_axes(
    [0.76, 0.08, 0.2, 0.2], projection="ipf", symmetry=sym
)
ax_ckey.plot_ipf_color_key(show_title=False)
ax_ckey.patch.set_facecolor("None")

directions = Vector3d(((1, 0, 0), (0, 1, 0), (0, 0, 1)))
n = directions.size

figsize = (4 * n * aspect_ratio, n * aspect_ratio)
fig, ax = plt.subplots(ncols=n, figsize=figsize)
for i, title in zip(range(n), ["X", "Y", "Z"]):
    ckey.direction = directions[i]
    rgb = ckey.orientation2color(xmap.rotations)
    ax[i].imshow(rgb.reshape(xmap.shape + (3,)))
    ax[i].set_title(f"IPF-{title}")
    ax[i].axis("off")
fig.subplots_adjust(wspace=0.02)

ref = ReciprocalLatticeVector(
    phase=xmap.phases[0], hkl=[[1, 1, 1], [2, 0, 0], [2, 2, 0], [3, 1, 1]]
)
ref = ref.symmetrise()
simulator = kp.simulations.KikuchiPatternSimulator(ref)
sim = simulator.on_detector(det, xmap.rotations.reshape(*xmap.shape))

# To remove existing markers
# del s.metadata.Markers

markers = sim.as_markers()
s.add_marker(markers, plot_marker=False, permanent=True)

maps_nav_rgb = kp.draw.get_rgb_navigator(rgb_x.reshape(xmap.shape + (3,)))

s.plot(maps_nav_rgb)
plt.show()

# Plot the phase map

xmap
xmap.plot()

fig = xmap.plot(remove_padding=True, return_figure=True)


plt.show()

Thank you in advance,
Andrea

[EffettiCollaterali]

Dear Håkon,
attached are three patterns extracted from the data of the previous comment. In particular, the first and the second ones are from austenite, the third is taken from a point of the grain that is indexed as Fe4N (by the Bruker software).
Also, I checked for the PCs estimation. We have 0.481 -0.160 the distance from the detector is 16.42 mm the tilt angle is 7.7°, the working distance is 4.0 mm and the lamella tilt is -18°. I will try again with Kikuchypy using these starting PCs.

Best,

Andrea

Pattern.zip