Live notebook

You can run this notebook in a live session Binder 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
from orix.crystal_map import PhaseList

import kikuchipy as kp


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]:
Title: Pattern
SignalType: EBSD
Array Chunk
Bytes 549.32 MiB 54.93 MiB
Shape (50, 50|480, 480) (25,10|480,480)
Count 11 Tasks 10 Chunks
Type uint8 numpy.ndarray

Navigation Axes

Signal Axes

50 50 480 480

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 | 105.73 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");
../_images/tutorials_pc_fit_plane_10_0.png

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]:
../_images/tutorials_pc_fit_plane_13_0.png

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:

  1. Extract an evenly spaced grid of patterns from the full dataset

  2. Get initial guesses of the PC for each grid pattern with Hough indexing (from PyEBSDIndex)

  3. Get initial guesses of the orientation for each grid pattern using the PCs with Hough indexing

  4. Refine the orientations and PCs simultaneously by matching experimental to dynamically simulated patterns (simulated with EMsoft)

  5. 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)
../_images/tutorials_pc_fit_plane_17_0.png

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 | 107.38 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,
)
../_images/tutorials_pc_fit_plane_21_0.png

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]:
phase = mp.phase

lat = phase.structure.lattice
lat.setLatPar(lat.a * 10, lat.b * 10, lat.c * 10)

print(phase)
print(phase.structure)
<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 | 106.01 ms
[########################################] | 100% Completed | 106.02 ms

Get an initial EBSD detector to store PC values in

[14]:
det_grid = s_grid.detector.deepcopy()

print(det_grid)
print(det_grid.sample_tilt)
EBSDDetector(shape=(480, 480), pc=(0.5, 0.5, 0.5), sample_tilt=70.0, tilt=0.0, azimuthal=0.0, binning=1.0, px_size=1.0 um)
70.0

We need an EBSDIndexer to use PyEBSDIndex. We can obtain an indexer by passing a PhaseList to EBSDDetector.get_indexer()

[15]:
pl = PhaseList(phase)
pl
[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(pl, 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]:
det_grid = s_grid.hough_indexing_optimize_pc(
    pc0=[0.52, 0.16, 0.49],  # Initial guess based on previous experiments
    indexer=indexer,
    batch=True,
    method="PSO",
    search_limit=0.1,
)

print(det_grid.pc_flattened.mean(axis=0))
print(det_grid.pc_flattened.std(0))

PC found: [**********] 25/25  global best:0.294  PC opt:[0.489  0.1323 0.4971]
[0.51985927 0.15997982 0.49042375]
[0.02294627 0.01848416 0.00770837]

Plot the PCs

[18]:
det_grid.plot_pc("map")
det_grid.plot_pc("scatter", annotate=True)
../_images/tutorials_pc_fit_plane_37_0.png
../_images/tutorials_pc_fit_plane_37_1.png

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(pl, nBands=6)
[20]:
xmap_grid = s_grid.hough_indexing(pl, indexer=indexer, verbose=2)

print(xmap_grid)
print(xmap_grid.fit.mean())
Hough indexing with PyEBSDIndex information:
  PyOpenCL: True
  Projection center (Bruker, mean): (0.5199, 0.16, 0.4904)
  Indexing 25 pattern(s) in 1 chunk(s)
Radon Time: 0.035934791987529024
Convolution Time: 0.02500349999172613
Peak ID Time: 0.026477624982362613
Band Label Time: 0.02196350001031533
Total Band Find Time: 0.11099612500402145
Band Vote Time:  0.00416962499730289
  Indexing speed: 197.48905 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.22957438
../_images/tutorials_pc_fit_plane_40_1.png

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)
../_images/tutorials_pc_fit_plane_42_0.png

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()
../_images/tutorials_pc_fit_plane_45_0.png

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]:
xmap_grid_ref, det_grid_ref = s_grid.refine_orientation_projection_center(
    xmap=xmap_grid,
    detector=det_grid,
    master_pattern=mp,
    energy=20,
    signal_mask=signal_mask,
    method="LN_NELDERMEAD",
    trust_region=[5, 5, 5, 0.05, 0.05, 0.05],
    rtol=1e-5,
    chunk_kwargs=dict(chunk_shape=1),
)
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 | 17.07 ss
Refinement speed: 1.46405 patterns/s

Check the average NCC score, number of evaluations, and PC and the standard deviation of the average PC

[24]:
print(xmap_grid_ref.scores.mean())
print(xmap_grid_ref.num_evals.mean())
print(det_grid_ref.pc_average)
print(det_grid_ref.pc_flattened.std(0))
0.3331923633813858
263.92
[0.5191021  0.15505618 0.48674104]
[0.0236548  0.02038152 0.0072576 ]

Check the distribution of refined PCs

[25]:
det_grid_ref.plot_pc("map")
det_grid_ref.plot_pc("scatter", annotate=True)
../_images/tutorials_pc_fit_plane_51_0.png
../_images/tutorials_pc_fit_plane_51_1.png

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 = det_grid_ref.estimate_xtilt(detect_outliers=False, degrees=True)
print(xtilt)
20.20873583715034
../_images/tutorials_pc_fit_plane_53_1.png

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(shape=(480, 480), pc=(0.52, 0.156, 0.486), sample_tilt=70.583, tilt=0.0, azimuthal=0.0, binning=1.0, px_size=1.0 um)
70.58272465784889
[0.0014476  0.00177993 0.00104693]
../_images/tutorials_pc_fit_plane_59_1.png
[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(shape=(480, 480), pc=(0.52, 0.156, 0.486), sample_tilt=70.855, tilt=0.0, azimuthal=0.0, binning=1.0, px_size=1.0 um)
70.85500334159653
[0.00130097 0.00143733 0.00106924]
../_images/tutorials_pc_fit_plane_60_1.png

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]:
xmap_grid_ref2 = s_grid.refine_orientation(
    xmap=xmap_grid_ref,
    detector=det_ref_proj_grid,
    master_pattern=mp,
    energy=20,
    signal_mask=signal_mask,
    method="LN_NELDERMEAD",
    trust_region=[5, 5, 5],
    rtol=1e-5,
    chunk_kwargs={"chunk_shape": 1},
)

print(xmap_grid_ref2.scores.mean())
Refinement information:
  Method: LN_NELDERMEAD (local) from NLopt
  Trust region (+/-): [5 5 5]
  Relative tolerance: 1e-05
Refining 25 orientation(s):
[########################################] | 100% Completed | 3.53 sms
Refinement speed: 7.06702 patterns/s
0.33362284004688264

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]:
g = ReciprocalLatticeVector.from_min_dspacing(phase.deepcopy())

# Ensure a complete unit cell (potentially changes the number of atoms)
g.sanitise_phase()

g.calculate_structure_factor()
F = abs(g.structure_factor)
g = g[F > 0.4 * F.max()]

g.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

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 | 101.61 ms
Finding zone axes that are in some pattern:
[########################################] | 100% Completed | 106.53 ms
Calculating detector coordinates for bands and zone axes:
[########################################] | 100% Completed | 106.46 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 | 105.76 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)
../_images/tutorials_pc_fit_plane_73_0.png

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 central 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.3467962169190736
../_images/tutorials_pc_fit_plane_76_1.png

We see that it is the least misoriented to the closest patterns, while the misorientation increases outwards radially. The first column of patterns are apparently all misorientation about 0.3\(^{\circ}\) to the central pattern. All misorientation angles are < 1\(^{\circ}\).