Indexing austenite and ferrite

[EffettiCollaterali]

Dear all I'm trying to index an austenite and ferrite sample.
I have a tiny map that I'm using for the tests. Here the data in HDF5 format.
https://drive.google.com/file/d/1tpZKcrfUBRTO-azx2gedmvJ1S01rpN8D/view?usp=sharing
I expect a central austenite grain and the outer part to be ferrite.
However, I am not able to get reasonable results. I tried to increase the number of reflectors decreasing the threshold on the structure factor, and increasing the number of bands. Also, I tried to downsize the pattern and to apply or not the background correction.
Any advice to improve the indexing would be appreciated.
Best,
Andrea

Phase map from bruker (the test is the map in the rectangular part (green austenite, red ferrite))
PhaseMap.pdf

Current input

`import matplotlib.pyplot as plt
import numpy as np
import pyvista
import skimage.exposure as ske
import skimage.transform as skt
import hyperspy
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
print('Load the data')

s=kp.load("New_NR_240_smallest.hdf5")
#s = kp.data.nickel_ebsd_large(allow_download=True)
#print(s.static_background)

#Plot the data
print('Plot the data')

s.plot()
plt.savefig('data.png')
plt.show()

s.plot()

#Remove the static background
print('Remove the static background')

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

print('Finished')

#Remove the dynamic background
print('Remove the dynamic background')

s.remove_dynamic_background()
s.plot()
plt.show()
print('Finished')

#Downsize the patterns
print('Downsize the patterns')

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

#s = s2
print('Finished')

#Print image quality
print('Print image quality')
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")

print('Finished')

#Set the detector
print('Set the detector')

det = kp.detectors.EBSDDetector(
shape=s.axes_manager.signal_shape[::-1],
pc=[0.5, 0.16, 0.5],
px_size=197.5, # Microns
binning=1,
tilt=6.924999999,
sample_tilt=69.99999,
)

det.pc_tsl()

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

#Plot the Virtual BSE
print('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.savefig('Virtual_BSE.png')
plt.show()
print('Finished')

#Plot the image quality
print('Plot the image quality')

maps_iq = s.get_image_quality()
print(maps_iq.mean())

plt.imshow(maps_iq, cmap="gray")
plt.tight_layout()
print('Finished')

#Plot the dot product
print('Plot the dot product')

adp = s.get_average_neighbour_dot_product_map()

s.xmap.plot(
adp.ravel(),
cmap="gray",
colorbar=True,
colorbar_label="Average dot product",
remove_padding=True,
)
print('Finished')

Extract data, also getting the grid positions

print('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,
)
print('Finished')

print('Make the phase list')

phase_list = PhaseList(
names=["Ferrite", "Austenite"],
space_groups=[229, 225],
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)
),
]
)

print('Make the reflectors')

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 = rlv[rlv.allowed]
rlv.sanitise_phase()
rlv = rlv.unique(use_symmetry=True).symmetrise()
rlv.calculate_structure_factor("lobato")
F = abs(rlv.structure_factor)
#print(F)
rlv = rlv[F > 0.4 * F.max()]
rlv.print_table()
reflectors.append(rlv)
print('Finished')

print('Sent the indexer')

indexer = det.get_indexer(
phase_list=phase_list,
reflectors=reflectors,
nBands = 10,
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)
#print(indexer.phaselist[2].latticeparameter)
#print(indexer.phaselist[2].polefamilies)
print('Finished')

print('Find the PCs')

det = s_grid.hough_indexing_optimize_pc(
pc0=[0.517, 0.169, 0.5], # Initial guess based on previous experiments
indexer=indexer,
batch=True,
method="PSO",
search_limit=0.02,
)
print('Finished')

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 = 12,
tSigma = 2.0,
rSigma = 2.0
)
indexer.PC

xmap = s.hough_indexing(phase_list=phase_list, chunksize=20, 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))

Plot overlayed images

fig = xmap.plot(rgb, overlay=maps_iq.ravel(), return_figure=True)

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.savefig('RGB.png')
plt.show()

Plot the phase map

xmap

xmap.plot()

fig = xmap.plot(remove_padding=True, return_figure=True)
plt.savefig('Phases.png')
plt.show()

`