This notebook is part of the kikuchipy documentation https://kikuchipy.org. Links to the documentation won't work from the notebook.

Extrapolate projection centers from a mean¶

In this tutorial, we will extrapolate a plane of projection centers (PCs) from a mean PC, as determined from patterns spread out across the sample region of interest. This is an alternative to fitting a plane to PCs, as done in the tutorial Fit a plane to selected projection centers. As a verification of the extrapolated PCs, we will compare them to the PCs obtained from fitting a plane to the PCs as is done in that tutorial.

We'll start by importing the necessary libraries

In [1]:

%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np

from diffsims.crystallography import ReciprocalLatticeVector
import hyperspy.api as hs
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¶

We will use nine calibration patterns from recrystallized nickel. The patterns have the full (480, 480) px$^2$ resolution of the NORDIF UF-1100 detector they are acquired on. These patterns are always acquired from spread out sample positions across the region of interest (ROI) in order to calibrate the PCs prior to indexing the full dataset.

Read the calibration patterns

In [2]:

s_cal = kp.data.ni_gain_calibration(1, allow_download=True)
s_cal

Out[2]:

<EBSD, title: Calibration patterns, dimensions: (9|480, 480)>

Get information read with NORDIF calibration patterns

In [3]:

omd = s_cal.original_metadata.as_dictionary()

Plot coordinates of calibration patterns on the secondary electron area overview image (part of the dataset), highlighting the region of interest (ROI)

In [4]:

kp.draw.plot_pattern_positions_in_map(
    rc=omd["calibration_patterns"]["indices_scaled"],
    roi_shape=omd["roi"]["shape_scaled"],
    roi_origin=omd["roi"]["origin_scaled"],
    area_shape=omd["area"]["shape_scaled"],
    area_image=omd["area_image"],
    color="w",
)

Plot the nine calibration patterns

In [5]:

fig = plt.figure(figsize=(10, 10))
_ = hs.plot.plot_images(
    s_cal,
    per_row=3,
    axes_decor=None,
    colorbar=False,
    label=None,
    fig=fig,
)
for i, ax in enumerate(fig.axes):
    ax.axis("off")
    ax.text(5, 10, str(i), va="top", ha="left", c="w")
fig.subplots_adjust(wspace=0.01, hspace=0.01)

We know that all patterns are of nickel. To get a description of nickel, we could create a Phase manually. However, we will later on use a dynamically simulated EBSD master pattern of nickel (created with EMsoft), which is loaded with a Phase. We will use this in the remaining analysis.

In [6]:

mp = kp.data.ebsd_master_pattern(
    "ni", allow_download=True, projection="lambert", energy=20
)
mp

Out[6]:

<EBSDMasterPattern, title: ni_mc_mp_20kv, dimensions: (|1001, 1001)>

Extract the phase, and change the lattice parameter unit from nm to Ångström

In [7]:

phase = mp.phase

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

In [8]:

print(phase)
print(phase.structure)

<name: ni. space group: Fm-3m. point group: m-3m. proper point group: 432. color: tab:blue>
lattice=Lattice(a=3.5236, b=3.5236, c=3.5236, alpha=90, beta=90, gamma=90)
28   0.000000 0.000000 0.000000 1.0000

Estimate PCs with Hough indexing¶

We will estimate PCs for the nine calibration 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, and can be installed with both pip and conda (from conda-forge). To install PyEBSDIndex, see their installation instructions.

PyEBSDIndex supports indexing face centered and body centered cubic (FCC and BCC) materials.

Improve signal-to-noise ratio by removing the static and dynamic background

In [9]:

s_cal.remove_static_background()
s_cal.remove_dynamic_background()

[########################################] | 100% Completed | 101.17 ms
[########################################] | 100% Completed | 101.52 ms

We need an EBSDIndexer to use PyEBSDIndex. We can obtain an indexer by passing a PhaseList to EBSDDetector.get_indexer(). Therefore, we need an initial EBSD detector

In [10]:

det_cal = s_cal.detector

print(det_cal)
print(det_cal.sample_tilt)

EBSDDetector (480, 480), px_size 1 um, binning 1, tilt 0.0, azimuthal 0.0, pc (0.5, 0.5, 0.5)
70.0

In [11]:

phase_list = PhaseList(phase)
phase_list

Out[11]:

Id  Name  Space group  Point group  Proper point group     Color
 0    ni        Fm-3m         m-3m                 432  tab:blue

In [12]:

indexer = det_cal.get_indexer(phase_list)
print(indexer.phaselist)

['FCC']

We estimate the PC of each pattern with Hough indexing, and plot both the mean and standard deviation of the resulting PCs. (We will "overwrite" the existing detector variable.)

In [13]:

det_cal = s_cal.hough_indexing_optimize_pc(
    pc0=[0.4, 0.2, 0.5],  # Initial guess based on previous experiments
    indexer=indexer,
    batch=True,
)

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

[0.42077662 0.2182365  0.49993497]
[0.0072175  0.00801777 0.00784516]

Plot the PCs

In [14]:

det_cal.plot_pc("scatter", annotate=True)

Unfortunately, we do not recognize the spatial distribution from the overview image above, and the expected inverse relation between (PCz, PCy) is not present... We need better indexing of the patterns, which we get by refining the PCs using pattern matching with dynamical simulations. These simulations are created with EMsoft.

First, we need an initial guess for the orientations, using Hough indexing via EBSD.hough_indexing(). We will use the mean PC for all patterns

In [15]:

indexer.PC = det_cal.pc_average

In [16]:

xmap_hi = s_cal.hough_indexing(
    phase_list=phase_list, indexer=indexer, verbose=2
)

Hough indexing with PyEBSDIndex information:
  PyOpenCL: True
  Projection center: (0.4208, 0.2182, 0.4999)
  Indexing 9 pattern(s) in 1 chunk(s)
Radon Time: 0.03481806100171525
Convolution Time: 0.004420306999236345
Peak ID Time: 0.0024418800021521747
Band Label Time: 0.04934866500116186
Total Band Find Time: 0.09106226299991249

Band Vote Time:  0.020639727998059243
  Indexing speed: 30.69000 patterns/s

Refine PCs with pattern matching¶

Refine the PCs (and orientations) using the Nelder-Mead implementation from NLopt

In [17]:

xmap_ref, det_ref = s_cal.refine_orientation_projection_center(
    xmap=xmap_hi,
    detector=det_cal,
    master_pattern=mp,
    energy=20,
    method="LN_NELDERMEAD",
    trust_region=[5, 5, 5, 0.05, 0.05, 0.05],
    rtol=1e-7,
    # A pattern per iteration to use all CPUs
    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-07
Refining 9 orientation(s) and projection center(s):
[########################################] | 100% Completed | 79.66 ss
Refinement speed: 0.11297 patterns/s

Inspect some refinement statistics

In [18]:

print("Score mean:                ", xmap_ref.scores.mean())
print("Mean number of evaluations:", xmap_ref.num_evals.mean())
print("PC mean:     ", det_ref.pc.mean(0))
print("PC std:      ", det_ref.pc.std(0))
print("PC max. diff:", abs(det_cal.pc - det_ref.pc).max(0))

ori_diff = xmap_hi.orientations.angle_with(
    xmap_ref.orientations, degrees=True
).max()
print("Ori. max. diff [deg]:", ori_diff)

Score mean:                 0.5265933639473386
Mean number of evaluations: 536.7777777777778
PC mean:      [0.41998311 0.21323872 0.50126773]
PC std:       [0.00269245 0.00188403 0.00060997]
PC max. diff: [0.00819502 0.01585595 0.01511784]
Ori. max. diff [deg]: 1.401081894104408

In [19]:

det_ref.plot_pc("scatter", annotate=True)

The diagonals 5-1-0-3-7 and 8-4-0-2-6 seen in the overview image above should be replicated in the (PCx, PCy) scatter plot. We see that 5-1-0-3-7 and 8-0-6 align as expected, but the PC values from the 2nd and 4th patterns do not lie in the expected range. Fitting a plane to all these values might not work to our satisfaction, so we will exclude the 2nd and 4th PC values when fitting a plane to the seven remaining PC values. The plane will have PC values for all points in the ROI in the overview image above.

In [20]:

is_outlier = np.zeros(det_ref.navigation_size, dtype=bool)
is_outlier[[2, 4]] = True

Get PC plane by fitting¶

The fitting is done by finding a transformation function which takes 2D sample coordinates and gives PC values for those coordinates. Both an affine and a projective transformation function is supported, following Winkelmann et al. (2020). By passing 2D indices of all ROI map points and of the points where the nine calibration patterns were obtained, EBSDDetector.fit_pc() returns a new detector with PC values for all map points. We will use the maximum difference between the above refined PC values and the corresponding fitted PC values as a measure of how good the fitted PC values are.

In [21]:

pc_indices = omd["calibration_patterns"]["indices_scaled"]
pc_indices -= omd["roi"]["origin_scaled"]
pc_indices = pc_indices.T

map_indices = np.indices(omd["roi"]["shape_scaled"])
print("Full map shape (n rows, n columns):", omd["roi"]["shape_scaled"])

Full map shape (n rows, n columns): (149, 200)

Fit PC values using the affine transformation function

In [22]:

det_fit_aff = det_ref.fit_pc(
    pc_indices=pc_indices,
    map_indices=map_indices,
    transformation="affine",
    is_outlier=is_outlier,
)

print(det_fit_aff.pc_average)
print(det_fit_aff.sample_tilt)

[0.41979185 0.21355389 0.50146853]
72.41196011494577

Fit PC values using the projective transformation function

In [23]:

det_fit_proj = det_ref.fit_pc(
    pc_indices=pc_indices,
    map_indices=map_indices,
    transformation="projective",
    is_outlier=is_outlier,
)

print(det_fit_proj.pc_average)
print(det_fit_proj.sample_tilt)

[0.4197901  0.21355393 0.50146771]
73.03021604645001

Compare PC differences of all but the PC values from the 2nd and 4th patterns

In [24]:

pc_indices2 = pc_indices.T[~is_outlier].T

# Refined PC values as a reference (ground truth)
pc_ref = det_ref.pc[~is_outlier]

# Difference in PC values from the affine transformation function
pc_diff_aff = det_fit_aff.pc[tuple(pc_indices2)] - pc_ref
pc_diff_aff_max = abs(pc_diff_aff).max(axis=0)
print(pc_diff_aff_max)

# Difference in PC values from the projective transformation function
pc_diff_proj = det_fit_proj.pc[tuple(pc_indices2)] - pc_ref
pc_diff_proj_max = abs(pc_diff_proj).max(axis=0)
print(pc_diff_proj_max)

# Is the difference from the affine transformation function greater than
# that from the projective one for (PCx, PCy, PCz)?
print(pc_diff_aff_max > pc_diff_proj_max)

[0.00063385 0.00036217 0.0002709 ]
[0.00067374 0.00033036 0.00029519]
[False  True False]

Get PC plane by extrapolation¶

Instead of fitting a plane to several PCs, we can extrapolate an average PC using EBSDDetector.extrapolate_pc(). To do this we need to know the detector pixel size and map step sizes, both given in the same unit.

In [25]:

det_ext = det_ref.extrapolate_pc(
    pc_indices=pc_indices,
    navigation_shape=omd["roi"]["shape_scaled"],
    step_sizes=(1.5, 1.5),  # um
    shape=det_cal.shape,
    px_size=70,  # In um. This is unique for every detector model!
    is_outlier=is_outlier,
)

print(det_ext.pc_average)
print(det_ext.sample_tilt)

[0.41887152 0.21397264 0.50132107]
70.0

Difference in PC values between the refined PCs and the corresponding extrapolated PCs

In [26]:

pc_diff_ext = det_ext.pc[tuple(pc_indices2)] - pc_ref
pc_diff_ext_max = abs(pc_diff_ext).max(axis=0)
print(pc_diff_ext_max)
print(pc_diff_ext_max > pc_diff_aff_max)

[0.00163081 0.00080355 0.0005684 ]
[ True  True  True]

The extrapolated PCs deviate more from the refined PCs than the PCs obtained from fitting, although the difference is small.

In [27]:

det_ext.plot_pc()

Validate extrapolated PCs¶

As a final check of the difference in PCs, we can plot the geometrical simulations on top of the patterns using the refined orientation but the three different PCs.

In [28]:

rlv = ReciprocalLatticeVector.from_min_dspacing(phase.deepcopy())

rlv.sanitise_phase()  # "Fill atoms in unit cell"

rlv.calculate_structure_factor()

structure_factor = abs(rlv.structure_factor)
rlv = rlv[structure_factor > 0.5 * structure_factor.max()]

rlv.print_table()

 h k l      d     |F|_hkl   |F|^2   |F|^2_rel   Mult 
 1 1 1    2.034    11.8     140.0     100.0      8   
 2 0 0    1.762    10.4     108.2      77.3      6   
 2 2 0    1.246     7.4     55.0       39.3      12  
 3 1 1    1.062     6.2     38.6       27.6      24

In [29]:

simulator = kp.simulations.KikuchiPatternSimulator(rlv)

We will use a signal mask during refinement to exclude pattern intensities in the detector corners

In [30]:

signal_mask = ~kp.filters.Window("circular", det_cal.shape).astype(bool)

Using PCs from the affine transformation function

In [31]:

det_aff_cal = det_fit_aff.deepcopy()
det_aff_cal.pc = det_aff_cal.pc[tuple(pc_indices)]

In [32]:

xmap_aff = s_cal.refine_orientation(
    xmap=xmap_hi,
    detector=det_aff_cal,
    master_pattern=mp,
    energy=20,
    signal_mask=signal_mask,
    method="LN_NELDERMEAD",
    trust_region=[5, 5, 5],
    rtol=1e-7,
    # A pattern per iteration to use all CPUs
    chunk_kwargs=dict(chunk_shape=1),
)

Refinement information:
  Method: LN_NELDERMEAD (local) from NLopt
  Trust region (+/-): [5 5 5]
  Relative tolerance: 1e-07
Refining 9 orientation(s):
[########################################] | 100% Completed | 13.64 ss
Refinement speed: 0.65973 patterns/s

In [33]:

sim_aff = simulator.on_detector(det_aff_cal, xmap_aff.rotations)

Finding bands that are in some pattern:
[########################################] | 100% Completed | 101.26 ms
Finding zone axes that are in some pattern:
[########################################] | 100% Completed | 101.34 ms
Calculating detector coordinates for bands and zone axes:
[########################################] | 100% Completed | 101.74 ms

Using PCs from the projective transformation function

In [34]:

det_proj_cal = det_fit_proj.deepcopy()
det_proj_cal.pc = det_proj_cal.pc[tuple(pc_indices)]

In [35]:

xmap_proj = s_cal.refine_orientation(
    xmap=xmap_hi,
    detector=det_proj_cal,
    master_pattern=mp,
    energy=20,
    signal_mask=signal_mask,
    method="LN_NELDERMEAD",
    trust_region=[5, 5, 5],
    rtol=1e-7,
    # A pattern per iteration to use all CPUs
    chunk_kwargs=dict(chunk_shape=1),
)

Refinement information:
  Method: LN_NELDERMEAD (local) from NLopt
  Trust region (+/-): [5 5 5]
  Relative tolerance: 1e-07
Refining 9 orientation(s):
[########################################] | 100% Completed | 14.40 ss
Refinement speed: 0.62480 patterns/s

In [36]:

sim_proj = simulator.on_detector(det_proj_cal, xmap_proj.rotations)

Finding bands that are in some pattern:
[########################################] | 100% Completed | 101.55 ms
Finding zone axes that are in some pattern:
[########################################] | 100% Completed | 101.89 ms
Calculating detector coordinates for bands and zone axes:
[########################################] | 100% Completed | 101.80 ms

Using the extrapolated PCs

In [37]:

det_ext_cal = det_ext.deepcopy()
det_ext_cal.pc = det_ext_cal.pc[tuple(pc_indices)]

In [38]:

xmap_ext = s_cal.refine_orientation(
    xmap=xmap_hi,
    detector=det_ext_cal,
    master_pattern=mp,
    energy=20,
    signal_mask=signal_mask,
    method="LN_NELDERMEAD",
    trust_region=[5, 5, 5],
    rtol=1e-7,
    # A pattern per iteration to use all CPUs
    chunk_kwargs=dict(chunk_shape=1),
)

Refinement information:
  Method: LN_NELDERMEAD (local) from NLopt
  Trust region (+/-): [5 5 5]
  Relative tolerance: 1e-07
Refining 9 orientation(s):
[########################################] | 100% Completed | 13.73 ss
Refinement speed: 0.65522 patterns/s

In [39]:

sim_ext = simulator.on_detector(det_ext_cal, xmap_ext.rotations)

Finding bands that are in some pattern:
[########################################] | 100% Completed | 101.97 ms
Finding zone axes that are in some pattern:
[########################################] | 100% Completed | 101.36 ms
Calculating detector coordinates for bands and zone axes:
[########################################] | 100% Completed | 102.20 ms

Compare normalized cross-correlation scores and number of evaluations (iterations)

In [40]:

print("Scores\n------")
print(f"Affine:       {xmap_aff.scores.mean():.7f}")
print(f"Projective:   {xmap_proj.scores.mean():.7f}")
print(f"Extrapolated: {xmap_ext.scores.mean():.7f}\n")

print("Number of evaluations\n---------------------")
print(f"Affine:       {xmap_aff.num_evals.mean():.1f}")
print(f"Projective:   {xmap_proj.num_evals.mean():.1f}")
print(f"Extrapolated: {xmap_ext.num_evals.mean():.1f}")

Scores
------
Affine:       0.5692966
Projective:   0.5692996
Extrapolated: 0.5685377

Number of evaluations
---------------------
Affine:       168.4
Projective:   176.7
Extrapolated: 169.3

Plot Kikuchi bands on top of patterns for the solutions using the affine transformed PCs (red), projective transformed PCs (blue) and extrapolated PCs (white)

In [41]:

fig, axes = plt.subplots(ncols=3, nrows=3, figsize=(12, 12))
for i, ax in enumerate(axes.ravel()):
    ax.imshow(s_cal.inav[i].data, cmap="gray")
    ax.axis("off")

    # Affine
    lines = sim_aff.as_collections(
        i, lines_kwargs=dict(linewidth=4, alpha=0.4)
    )[0]
    ax.add_collection(lines)

    # Projective
    lines = sim_proj.as_collections(
        i, lines_kwargs=dict(color="b", linewidth=3, alpha=0.4)
    )[0]
    ax.add_collection(lines)

    # Extrapolated
    lines = sim_ext.as_collections(
        i, lines_kwargs=dict(color="w", linewidth=1, alpha=0.4)
    )[0]
    ax.add_collection(lines)

    ax.text(5, 10, i, c="w", va="top", ha="left", fontsize=20)
fig.subplots_adjust(wspace=0.01, hspace=0.01)

As expected from the intermediate results above (similar PCs and mean NCC score), all PCs produce visually identical geometrical simulations. However, the orientations are still slightly different, as we see by calculating the misorientation angles between the orientations obtained using the extrapolated PCs and the fitted PCs

In [42]:

xmap_proj.orientations.angle_with(xmap_ext.orientations, degrees=True)

Out[42]:

array([3.00016517, 2.99419537, 2.99722527, 3.00370345, 3.00205733,
       2.99066348, 2.99660149, 3.00870014, 3.01203182])

The misorientation of 3$^{\circ}$ most likely comes from the difference in applied sample tilt between the detector used for the projected and extrapolated PCs. Since these are experimental data, it's difficult to say which sample tilt is more correct, although the nominal sample tilt from the microscope was 70$^{\circ}$.