Live notebook
You can run this notebook in a live session 
or view it on Github.
Fit a plane to selected projection centers#
In this tutorial, we will fit a plane to projection centers (PCs). The PCs are determined from a (5, 5) grid of EBSD patterns from a single crystal silicon wafer. The resulting fitted plane contains one PC per pattern in the dataset from which the grid is obtained. This plane of PCs can be used in subsequent indexing of the dataset. To relate the sample positions of the grid to the PCs, we will test both an affine transformation and a projective transformation. These transformations are based on the work by [Winkelmann et al., 2020].
We’ll start by importing necessary libraries
[1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import skimage.filters as skf
from diffsims.crystallography import ReciprocalLatticeVector
import kikuchipy as kp
from orix.crystal_map import PhaseList
plt.rcParams.update(
{
"figure.facecolor": "w",
"figure.dpi": 75,
"figure.figsize": (8, 8),
"font.size": 15,
}
)
Load and inspect data#
Load data lazily (~500 MB) and inspect its shape and step size
[2]:
s = kp.data.si_wafer(allow_download=True, lazy=True)
s
[2]:
|
|
We see that we have (50, 50) patterns of shape (480, 480). This is the full pattern resolution of the NORDIF UF-420 detector the patterns are acquired on. In the axes manager…
[3]:
print(s.axes_manager)
<Axes manager, axes: (50, 50|480, 480)>
Name | size | index | offset | scale | units
================ | ====== | ====== | ======= | ======= | ======
x | 50 | 0 | 0 | 40 | um
y | 50 | 0 | 0 | 40 | um
---------------- | ------ | ------ | ------- | ------- | ------
dx | 480 | 0 | 0 | 1 | um
dy | 480 | 0 | 0 | 1 | um
… we see that the nominal step sizes are \((\Delta x, \Delta y) = (40, 40)\) \(\mu\)m. The scan region of interest (ROI) therefore covers an area of \((2 \times 2)\) \(\mu\)m\(^2\). The NORDIF UF-420 (unbinned) detector pixel size is about 90 \(\mu\)m. This means that we expect the PC to shift about 2000 / 90 \(\approx\) 22 detector pixels when moving horizontally in the scan. To get a better understanding of the expected changes of the PC in the scan, we’ll look at the mean intensity map of the sample.
[4]:
s_mean_nav = s.mean(axis=(2, 3))
s_mean_nav.compute() # Need this call because data was loaded lazily
[########################################] | 100% Completed | 315.08 ms
Plot it, also annotating the orientation of the scan relative to the sample position in the chamber. We use the plot() method of the EBSD.xmap attribute for this.
[5]:
s.xmap.scan_unit = "um"
fig = s.xmap.plot(
s_mean_nav.data.ravel(),
colorbar=True,
colorbar_label="Mean intensity",
return_figure=True,
)
# Annotate
ax = fig.axes[0]
ax.text(25, 1, "Bottom", va="top", ha="center")
_ = ax.text(25, 49, "Top", va="bottom", ha="center");
The brighter and darker spots come from patterns acquired close to fiducial markers on the Si wafer. The vertical intensity gradient is a result of the intensity distribution on the detector shifting across the sample. Note the orientation of the ROI on the sample: to avoid depositing carbon ahead of the scanning beam, the sample is scanned from down the sample upwards, annotated on the map (“Bottom”, “Top”).
The Bruker PC convention is used internally in kikuchipy. The detector is viewed from the detector towards the sample. The convention measures (PCx, PCy, PCz) in the following manner:
PCx from left towards the right
PCy downwards from the top
PCz as the shortest distance from the beam-sample interaction position to the detector in fractions of the detector width.
See the reference frames tutorial for more details and conversions between conventions.
Given this definition, let’s add annotations of the expected changes to the PC (increasing in the direction of the arrow) to the map
[6]:
annotate_kw = {
"arrowprops": {"arrowstyle": "-|>", "lw": 2, "mutation_scale": 15}
}
xys = [(35, 25), (40, 35), (40, 15)]
xy_texts = [(49, 25), (40, 49), (40, 1)]
has = ["right", "center", "center"]
vas = ["center", "bottom", "top"]
labels = ["PCx", "PCy", "PCz"]
for xy, xytext, ha, va, label in zip(xys, xy_texts, has, vas, labels):
ax.annotate(label, xy=xy, xytext=xytext, ha=ha, va=va, **annotate_kw)
fig # Show figure again with new annotations
[6]:
Our goal is to fit a plane of (PCx, PCy, PCz) to this (50, 50) map with these variations. To this end, we perform the following actions:
Extract an evenly spaced grid of patterns from the full dataset
Get initial guesses of the PC for each grid pattern with Hough indexing (from PyEBSDIndex)
Get initial guesses of the orientation for each grid pattern using the PCs with Hough indexing
Refine the orientations and PCs simultaneously by matching experimental to dynamically simulated patterns (simulated with EMsoft)
Fit a plane to the refined PCs
Extract a grid of patterns#
Before we extract the grid of patterns and prepare these for indexing by background correction, we obtain a binary mask to remove parts of the pattern without Kikuchi diffraction (typically the corners). For this dataset, we get the mask via thresholding. The threshold value is found by using the minimum threshold algorithm implemented in scikit-image (note that they have other thresholding algorithms that might work better for other datasets).
[7]:
s_mean = s.mean((0, 1)) # (480, 480)
mean_data = s_mean.data.compute()
[8]:
# Threshold
threshold = skf.threshold_minimum(mean_data)
signal_mask = mean_data < threshold # Exclude pixels set to True
# Extract center pattern for plotting
p0 = s.inav[24, 24].data
fig, axes = plt.subplots(ncols=3, figsize=(12, 4))
for ax, arr, title in zip(
axes,
[p0, signal_mask, p0 * ~signal_mask],
["Pattern", "Mask", "Pattern * ~Mask"],
):
ax.imshow(arr, cmap="gray")
ax.axis("off")
ax.set(title=title)
fig.subplots_adjust(wspace=0)
The mask efficiently removes the parts of the pattern which contains no information of interest.
Extract the grid using EBSD.extract_grid()
[9]:
grid_shape = (5, 5)
s_grid, idx = s.extract_grid(grid_shape, return_indices=True)
nav_shape_grid = s_grid.axes_manager.navigation_shape[::-1]
sig_shape_grid = s_grid.axes_manager.signal_shape[::-1]
s_grid.compute()
s_grid
[########################################] | 100% Completed | 103.18 ms
[9]:
<EBSD, title: Pattern, dimensions: (5, 5|480, 480)>
Let’s plot the grid positions on the navigation map
[10]:
kp.draw.plot_pattern_positions_in_map(
rc=idx.reshape(2, -1).T,
roi_shape=(50, 50),
roi_image=s_mean_nav.data,
)
We know that all patterns are of silicon. To get a description of silicon, we could create a Phase manually. However, we will later on use a dynamically simulated EBSD master pattern of silicon (created with EMsoft), which is loaded with a Phase. We will use this in the remaining analysis.
[11]:
mp = kp.data.ebsd_master_pattern(
"si", allow_download=True, projection="lambert", energy=20
)
mp
[11]:
<EBSDMasterPattern, title: si_mc_mp_20kv, dimensions: (|1001, 1001)>
Extract the phase and change lattice constant unit from nm to Ångström
[12]:
<name: si. space group: Fd-3m. point group: m-3m. proper point group: 432. color: tab:blue>
lattice=Lattice(a=5.4307, b=5.4307, c=5.4307, alpha=90, beta=90, gamma=90)
14 0.000000 0.000000 0.000000 1.0000
Initial guess of PC using Hough indexing#
We will estimate PCs and do Hough indexing of the grid of patterns using PyEBSDIndex. See the Hough indexing tutorial for more details on the Hough indexing related commands below.
Note
PyEBSDIndex is an optional dependency of kikuchipy. It can be installed with pip or conda (from conda-forge). To install PyEBSDIndex, see their installation instructions.
Prepare the patterns for indexing by background removal and setting the masked out values to 0
[13]:
s_grid.remove_static_background()
s_grid.remove_dynamic_background()
s_grid = s_grid * ~signal_mask
[########################################] | 100% Completed | 102.00 ms
[########################################] | 100% Completed | 202.01 ms
Get an initial EBSD detector to store PC values in
EBSDDetector (480, 480), px_size 1 um, binning 1, tilt 0.0, azimuthal 0.0, pc (0.5, 0.5, 0.5)
70.0
We need an EBSDIndexer to use PyEBSDIndex. We can obtain an indexer by passing a PhaseList to EBSDDetector.get_indexer()
[15]:
phase_list = PhaseList(phase)
phase_list
[15]:
Id Name Space group Point group Proper point group Color
0 si Fd-3m m-3m 432 tab:blue
Get the indexer
[16]:
indexer = det_grid.get_indexer(phase_list, nBands=6)
We estimate the PC of each pattern with Hough indexing using EBSD.hough_indexing_optimize_pc(), and plot both the mean and standard deviation of the resulting PCs. (We will “overwrite” the existing detector variable in memory.)
[17]:
PC found: [**********] 25/25 global best:0.274 PC opt:[0.4899 0.1319 0.4973]
[0.52006439 0.1596286 0.4907931 ]
[0.02247494 0.01867778 0.00789391]
Plot the PCs
[18]:
det_grid.plot_pc("map")
det_grid.plot_pc("scatter", annotate=True)
Looking at the maps, as expected:
PCx increases towards the left
PCy increases upwards
PCz increases downwards
Looking at the scatter plots:
We recognize the regular (5, 5) grid in the (PCx, PCy) plot.
There is an inverse relation between PCy and PCz: as PCy decreases, PCz increases. This is as expected. However, there is a significant spread in each row of points (5-9, 10-14, etc.).
These initial guesses seem OK and we could perhaps fit a plane to these values already. But, we can get a better fit by refining the PC values using pattern matching. To that end, we obtain initial guesses for the orientations via Hough indexing with EBSD.hough_indexing(). We get a new indexer with the optimized PCs and then index
[19]:
indexer = det_grid.get_indexer(phase_list, nBands=6)
[20]:
xmap_grid = s_grid.hough_indexing(
phase_list=phase_list, indexer=indexer, verbose=2
)
print(xmap_grid)
print(xmap_grid.fit.mean())
Hough indexing with PyEBSDIndex information:
PyOpenCL: True
Projection center (Bruker, mean): (0.5201, 0.1596, 0.4908)
Indexing 25 pattern(s) in 1 chunk(s)
Radon Time: 0.04952027699619066
Convolution Time: 0.005532572002266534
Peak ID Time: 0.0023458810028387234
Band Label Time: 0.06669590099772904
Total Band Find Time: 0.12412348898942582
Band Vote Time: 0.020521485988865606
Indexing speed: 147.07505 patterns/s
Phase Orientations Name Space group Point group Proper point group Color
0 25 (100.0%) si Fd-3m m-3m 432 tab:blue
Properties: fit, cm, pq, nmatch
Scan unit: um
0.20949268
Check the average pattern fit and confidence metric
[21]:
fig, axes = plt.subplots(ncols=2, figsize=(10, 3.5))
for ax, to_plot, label in zip(
axes, ["fit", "cm"], ["Pattern fit [deg]", "Confidence metric"]
):
im = ax.imshow(xmap_grid.get_map_data(to_plot))
fig.colorbar(im, ax=ax, label=label, pad=0.02)
ax.axis("off")
fig.subplots_adjust(wspace=0.1)
Most patterns have a low pattern fit and high confidence metric.
Let’s refine these initial guesses of orientations and PCs using the dynamically simulated master pattern.
Refine PCs using pattern matching#
Let’s look at the dynamical simulation (which is in the square Lambert projection, required for use in dictionary indexing and refinement in kikuchipy [and EMsoft])
[22]:
mp.plot()
Refine orientations and PCs using the Nelder-Mead optimization algorithm as implemented in NLopt. During optimization, orientations are represented as three Euler angles (\(\phi_1, \Phi, \phi_2\)). The angles are varied within a trust region in degrees of the orientation obtained from Hough indexing above. PCs are varied within a trust region in percent of the detector width (for PCx and PCz) or height (for PCy). Optimization of each pattern stops when the increase in normalized cross-correlation (NCC) score between iterations, obtained by comparing experimental to simulated pattern, is lower than \(10^{-5}\).
[23]:
Refinement information:
Method: LN_NELDERMEAD (local) from NLopt
Trust region (+/-): [5. 5. 5. 0.05 0.05 0.05]
Relative tolerance: 1e-05
Refining 25 orientation(s) and projection center(s):
[########################################] | 100% Completed | 114.76 s
Refinement speed: 0.21783 patterns/s
Check the average NCC score, number of evaluations, and PC and the standard deviation of the average PC
[24]:
0.33365421533584594
291.44
[0.51870718 0.15458022 0.48636834]
[0.02396429 0.02004163 0.00739694]
Check the distribution of refined PCs
[25]:
det_grid_ref.plot_pc("map")
det_grid_ref.plot_pc("scatter", annotate=True)
The grid is easily recognizable in all three scatter plots now. If there were any outliers, we could exclude them by passing a boolean array when fitting a plane. This array can either be created manually, or we can try to find outliers by robust fitting of a line to the (PCz, PCy) plot. Assuming the sample is perfectly tilted about the detector \(X_d\), the slope of the fitted line is the estimated tilt angle.
[26]:
xtilt, is_outlier = det_grid_ref.estimate_xtilt(
detect_outliers=True, degrees=True, return_outliers=True
)
print(xtilt)
20.550490605557012
The sample was tilted to 70\(^{\circ}\) from the horizontal, so we would expect an angle of about 20\(^{\circ}\) from the vertical.
We can now fit a plane to the remaining PCs.
Fit a plane to the refined PCs#
To fit a plane to the PCs, we must pass an array of all indices in the full map and the indices of the patterns which the PCs were estimated from. In the case of the grid patterns we’ve used here, the map indices were returned from the EBSD.extract_grid().
[27]:
nav_shape = s.axes_manager.navigation_shape[::-1]
map_indices = np.indices(nav_shape)
Using EBSDDetector.fit_pc(), we will fit a plane to the PCs using both an affine transformation and a projective transformation. The fit method automatically plots the three scatter plots above and a forth 3D plot. The experimental PC values are black circles and the fitted PCs are larger gray circles.
[28]:
det_ref_aff = det_grid_ref.fit_pc(
idx,
map_indices=map_indices,
transformation="affine",
)
print(det_ref_aff)
# Sample tilt
print(det_ref_aff.sample_tilt)
# Max. deviation between experimental and fitted PC
pc_diff_aff = det_grid_ref.pc - det_ref_aff.pc[tuple(idx)]
print(abs(pc_diff_aff.reshape(-1, 3)).mean(axis=0))
EBSDDetector (480, 480), px_size 1 um, binning 1, tilt 0.0, azimuthal 0, pc (0.52, 0.155, 0.486)
69.8010018609663
[0.00092969 0.0014497 0.00065325]
[29]:
det_ref_proj = det_grid_ref.fit_pc(
idx,
map_indices=map_indices,
transformation="projective",
)
print(det_ref_proj)
# Sample tilt
print(det_ref_proj.sample_tilt)
# Max. deviation between experimental and fitted PC
pc_diff_proj = det_grid_ref.pc - det_ref_proj.pc[tuple(idx)]
print(abs(pc_diff_proj.reshape(-1, 3)).mean(axis=0))
EBSDDetector (480, 480), px_size 1 um, binning 1, tilt 0.0, azimuthal 0, pc (0.52, 0.155, 0.486)
69.96587791952271
[0.00105687 0.00110206 0.00063361]
Both planes fit the experimental PCs well.
Validate fitted PCs#
Finally, we refine the (already refined) orientations of the grid patterns using the above fitted PCs. We validate our results by comparing geometrical simulations to the patterns.
[30]:
det_ref_proj_grid = det_ref_proj.deepcopy()
det_ref_proj_grid.pc = det_ref_proj_grid.pc[tuple(idx)]
[31]:
Refinement information:
Method: LN_NELDERMEAD (local) from NLopt
Trust region (+/-): [5 5 5]
Relative tolerance: 1e-05
Refining 25 orientation(s):
[########################################] | 100% Completed | 20.83 ss
Refinement speed: 1.19951 patterns/s
0.3322068232297897
To plot geometrical simulations on top of our experimental patterns, we create a KikuchiPatternSimulator using the silicon Phase from the master pattern above (see the geometrical EBSD simulations tutorial for more details)
[32]:
ref = ReciprocalLatticeVector.from_min_dspacing(phase)
# Ensure a complete unit cell (potentially changes the number of atoms)
ref.sanitise_phase()
ref.calculate_structure_factor()
F = abs(ref.structure_factor)
ref = ref[F > 0.4 * F.max()]
ref.print_table()
h k l d |F|_hkl |F|^2 |F|^2_rel Mult
1 1 1 3.135 18.4 338.2 100.0 8
2 2 0 1.920 15.0 224.0 66.2 12
4 0 0 1.358 9.2 83.7 24.8 6
3 1 1 1.637 8.5 72.4 21.4 24
[33]:
Get one simulation per pattern
[34]:
sim = simulator.on_detector(
det_ref_proj_grid,
xmap_grid_ref2.rotations.reshape(*xmap_grid_ref2.shape),
)
Finding bands that are in some pattern:
[########################################] | 100% Completed | 103.73 ms
Finding zone axes that are in some pattern:
[########################################] | 100% Completed | 101.46 ms
Calculating detector coordinates for bands and zone axes:
[########################################] | 100% Completed | 100.83 ms
Plot the geometrical simulations on top of the patterns (after normalization to a mean of 0 and standard deviation of 1)
[35]:
s_grid2 = s_grid.normalize_intensity(dtype_out="float32", inplace=False)
[########################################] | 100% Completed | 100.92 ms
[36]:
fig, axes = plt.subplots(*grid_shape, figsize=(15, 15))
for i in np.ndindex(grid_shape):
axes[i].imshow((s_grid2.data[i] * ~signal_mask), cmap="gray", vmin=-3, vmax=3)
axes[i].axis("off")
lines = sim.as_collections(i)[0]
axes[i].add_collection(lines)
idx1d = np.ravel_multi_index(i, grid_shape)
axes[i].text(5, 10, idx1d, c="w", va="top", ha="left", fontsize=20)
fig.subplots_adjust(wspace=0.01, hspace=0.01)
As expected, if we look at the first pattern (0), we see that the patterns shift upwards and to the left as we move down (see how the upper-most near-horizontal Kikuchi line disappears) and to the right (see how the left-most zone axis disappears) in the grid.
The patterns come from a single crystal. How misoriented is the estimated orientation for the center pattern to the other orientations?
[37]:
angles = xmap_grid_ref2.orientations[12].angle_with(
xmap_grid_ref2.orientations, degrees=True
)
print(angles.max())
xmap_grid_ref2.plot(angles, colorbar=True, colorbar_label="Misorientation angle [deg]")
0.7114465443203103
We see that it is the least misoriented to the closest patterns, while the misorientation increases outwards radially. All misorientation angles are < 1\(^{\circ}\).