This page was generated from doc/pattern_processing.ipynb.

Pattern processing

The raw EBSD signal can be empirically evaluated as a superposition of a Kikuchi diffraction pattern and a smooth background intensity. For pattern indexing, the latter intensity is usually undesirable, while for virtual backscatter electron (VBSE) imaging, this intensity can reveal topographical, compositional or diffraction contrast. This section details methods to enhance the Kikuchi diffraction pattern and manipulate detector intensities in patterns in an EBSD signal.

Most of the methods operating on EBSD objects use functions that operate on the individual patterns (numpy.ndarray). These single pattern functions are available in the kikuchipy.pattern module.

Let’s import the necessary libraries and read the Nickel EBSD test data set:

[1]:
# exchange inline for qt5 for interactive plotting from the pyqt package
%matplotlib inline

import matplotlib.pyplot as plt
plt.rcParams["font.size"] = 15
import numpy as np
import kikuchipy as kp


s = kp.data.nickel_ebsd_small()  # Use kp.load("data.h5") to load your own data
WARNING:hyperspy.api:The ipywidgets GUI elements are not available, probably because the hyperspy_gui_ipywidgets package is not installed.
WARNING:hyperspy.api:The traitsui GUI elements are not available, probably because the hyperspy_gui_traitsui package is not installed.

Most methods operate inplace (indicated in their docstrings), meaning they overwrite the patterns in the EBSD signal. If we instead want to keep the original signal and operate on a new signal, we can create a deepcopy() of the original signal. As an example here, we create a new EBSD signal from a small part of the original signal:

[2]:
s2 = s.deepcopy()
np.may_share_memory(s.data, s2.data)
[2]:
False

Background correction

Remove the static background

Effects which are constant, like hot pixels or dirt on the detector, can be removed by either subtracting or dividing by a static background via remove_static_background():

[3]:
s2.remove_static_background(operation="subtract", relative=True)

# _ means we don't want the output
_, ax = plt.subplots(figsize=(13, 6), ncols=2)
ax[0].imshow(s.inav[0, 0].data, cmap="gray")
ax[0].set_title("As acquired")
ax[1].imshow(s2.inav[0, 0].data, cmap="gray")
_ = ax[1].set_title("Static background removed")
Removing the static background:
[########################################] | 100% Completed |  0.7s
_images/pattern_processing_6_1.png

Here, the static background pattern is assumed to be stored as part of the signal metadata, which can be loaded via set_experimental_parameters(). The static background pattern can also be passed to the static_bg parameter. Passing relative=True (default) ensures that relative intensities between patterns are kept when they are rescaled after correction to fill the available data range. In this case, for a scan of data type uint8 with data range [0, 255], the highest pixel intensity in a scan is stretched to 255 (and the lowest to 0), while the rest is rescaled keeping relative intensities between patterns. With relative=False, all patterns are stretched to [0, 255].

The static background pattern intensities can be rescaled to each individual pattern’s intensity range before removal by passing scale_bg=True, which will result in the relative intensity between patterns to be lost (passing relative=True along with scale_bg=True is not allowed).

Remove the dynamic background

Uneven intensity in a static background subtracted pattern can be corrected by subtracting or dividing by a dynamic background obtained by Gaussian blurring. This so-called flat fielding is done with remove_dynamic_background(). A Gaussian window with a standard deviation set by std is used to blur each pattern individually (dynamic) either in the spatial or frequency domain, set by filter_domain. Blurring in the frequency domain is effectively accomplished by a low-pass Fast Fourier Transform (FFT) filter. The individual Gaussian blurred dynamic backgrounds are then subtracted or divided from the respective patterns, set by operation:

[4]:
s3 = s2.deepcopy()
s3.remove_dynamic_background(
    operation="subtract",  # Default
    filter_domain="frequency",  # Default
    std=8,  # Default is 1/8 of the pattern width
    truncate=4,  # Default
)

_, ax = plt.subplots(figsize=(13, 6), ncols=2)
ax[0].imshow(s2.inav[0, 0].data, cmap="gray")
ax[0].set_title("Static background removed")
ax[1].imshow(s3.inav[0, 0].data, cmap="gray")
_ = ax[1].set_title("Static + dynamic background removed")
Removing the dynamic background:
[########################################] | 100% Completed |  0.4s
_images/pattern_processing_9_1.png

The width of the Gaussian window is truncated at the truncated number of standard deviations. Output patterns are rescaled to fill the input patterns’ data type range.


Get the dynamic background

The Gaussian blurred pattern removed during dynamic background correction can be obtained as an EBSD signal by calling get_dynamic_background():

[5]:
bg = s.get_dynamic_background(filter_domain="frequency", std=8, truncate=4)
bg

_, ax = plt.subplots(figsize=(13, 6), ncols=2)
ax[0].imshow(s.inav[0, 0].data, cmap="gray")
ax[0].set_title("As acquired")
ax[1].imshow(bg.inav[0, 0].data, cmap="gray")
_ = ax[1].set_title("Dynamic background")
Getting the dynamic background:
[########################################] | 100% Completed |  0.1s
_images/pattern_processing_12_1.png

Average neighbour patterns

The signal-to-noise ratio in patterns in an EBSD signal s can be improved by averaging patterns with their closest neighbours within a window/kernel/mask with average_neighbour_patterns():

[6]:
s4 = s3.deepcopy()
s4.average_neighbour_patterns(window="gaussian", std=1)

_, ax = plt.subplots(figsize=(13, 6), ncols=2)
ax[0].imshow(s3.inav[0, 0].data, cmap="gray")
ax[0].set_title("Static + dynamic background removed")
ax[1].imshow(s4.inav[0, 0].data, cmap="gray")
_ = ax[1].set_title("After neighbour pattern averaging")
Averaging with the neighbour patterns:
[########################################] | 100% Completed |  0.1s
_images/pattern_processing_14_1.png

The array of averaged patterns \(g(n_{\mathrm{x}}, n_{\mathrm{y}})\) is obtained by spatially correlating a window \(w(s, t)\) with the array of patterns \(f(n_{\mathrm{x}}, n_{\mathrm{y}})\), here 4D, which is padded with zeros at the edges. As coordinates \(n_{\mathrm{x}}\) and \(n_{\mathrm{y}}\) are varied, the window origin moves from pattern to pattern, computing the sum of products of the window coefficients with the neighbour pattern intensities, defined by the window shape, followed by normalizing by the sum of the window coefficients. For a symmetrical window of shape \(m \times n\), this becomes [GW17]

\begin{equation} g(n_{\mathrm{x}}, n_{\mathrm{y}}) = \frac{\sum_{s=-a}^a\sum_{t=-b}^b{w(s, t) f(n_{\mathrm{x}} + s, n_{\mathrm{y}} + t)}} {\sum_{s=-a}^a\sum_{t=-b}^b{w(s, t)}}, \end{equation}

where \(a = (m - 1)/2\) and \(b = (n - 1)/2\). The window \(w\), a Window object, can be plotted:

[7]:
w = kp.filters.Window(window="gaussian", shape=(3, 3), std=1)
_, _, _ = w.plot(cmap="inferno")
_images/pattern_processing_16_0.png

Any 1D or 2D window with desired coefficients can be used. This custom window can be passed to the window parameter in average_neighbour_patterns() or Window as a numpy.ndarray or a dask.array.Array. Additionally, any window in scipy.signal.windows.get_window() passed as a string via window with the necessary parameters as keyword arguments (like std=1 for window="gaussian") can be used. To demonstrate the creation and use of an asymmetrical circular window (and the use of make_circular(), although we could create a circular window directly by calling window="circular" upon window initialization):

[8]:
w = kp.filters.Window(window="rectangular", shape=(5, 4))
w
[8]:
Window (5, 4) rectangular
[[1. 1. 1. 1.]
 [1. 1. 1. 1.]
 [1. 1. 1. 1.]
 [1. 1. 1. 1.]
 [1. 1. 1. 1.]]
[9]:
w.make_circular()
w
[9]:
Window (5, 4) circular
[[0. 0. 1. 0.]
 [0. 1. 1. 1.]
 [1. 1. 1. 1.]
 [0. 1. 1. 1.]
 [0. 0. 1. 0.]]
[10]:
s5 = s3.deepcopy()
s5.average_neighbour_patterns(w)
Averaging with the neighbour patterns:
[########################################] | 100% Completed |  0.1s
[11]:
fig, im, cbar = w.plot()
_images/pattern_processing_21_0.png

But this (5, 4) averaging window cannot be used with our (3, 3) navigation shape signal.

Note

Neighbour pattern averaging increases the virtual interaction volume of the electron beam with the sample, leading to a potential loss in spatial resolution. Averaging may in some cases, like on grain boundaries, mix two or more different diffraction patterns, which might be unwanted. See [WNL+15] for a discussion of this concern.


Adaptive histogram equalization

Enhancing the pattern contrast with adaptive histogram equalization has been found useful when comparing patterns for dictionary indexing [MDeGraefS+17]. With adaptive_histogram_equalization(), the intensities in the pattern histogram are spread to cover the available range, e.g. [0, 255] for patterns of uint8 data type:

[12]:
s6 = s3.deepcopy()
s6.adaptive_histogram_equalization(kernel_size=(15, 15))

_, ax = plt.subplots(figsize=(13, 6), ncols=2)
ax[0].imshow(s3.inav[0, 0].data, cmap="gray")
ax[0].set_title("Static + dynamic background removed")
ax[1].imshow(s6.inav[0, 0].data, cmap="gray")
_ = ax[1].set_title("After adaptive histogram equalization")
Adaptive histogram equalization:
[########################################] | 100% Completed |  0.1s
_images/pattern_processing_25_1.png

The kernel_size parameter determines the size of the contextual regions. See e.g. Fig. 5 in [JPDeGraef19], also available via EMsoft’s GitHub repository wiki, for the effect of varying kernel_size.


Filtering in the frequency domain

Filtering of patterns in the frequency domain can be done with fft_filter(). This method takes a spatial kernel defined in the spatial domain, or a transfer function defined in the frequency domain, in the transfer_function argument as a numpy.ndarray or a Window. Which domain the transfer function is defined in must be passed to the function_domain argument. Whether to shift zero-frequency components to the centre of the FFT can also be controlled via shift, but note that this is only used when function_domain="frequency".

Popular uses of filtering of EBSD patterns in the frequency domain include removing large scale variations across the detector with a Gaussian high pass filter, or removing high frequency noise with a Gaussian low pass filter. These particular functions are readily available via Window:

[13]:
pattern_shape = s.axes_manager.signal_shape[::-1]
w_low = kp.filters.Window(
    window="lowpass",
    cutoff=23,
    cutoff_width=10,
    shape=pattern_shape
)
w_high = kp.filters.Window(
    window="highpass",
    cutoff=3,
    cutoff_width=2,
    shape=pattern_shape
)

fig, ax = plt.subplots(figsize=(22, 6), ncols=3)
fig.subplots_adjust(wspace=0.05)
im0 = ax[0].imshow(w_low, cmap="gray")
ax[0].set_title("Lowpass filter", fontsize=22)
fig.colorbar(im0, ax=ax[0])
im1 = ax[1].imshow(w_high, cmap="gray")
ax[1].set_title("Highpass filter", fontsize=22)
fig.colorbar(im1, ax=ax[1])
im2 = ax[2].imshow(w_low * w_high, cmap="gray")
ax[2].set_title("Lowpass * highpass filter", fontsize=22)
_ = fig.colorbar(im2, ax=ax[2])
_images/pattern_processing_28_0.png

Then, to multiply the FFT of each pattern with this transfer function, and subsequently computing the inverse FFT (IFFT), we use fft_filter(), and remember to shift the zero-frequency components to the centre of the FFT:

[14]:
s7 = s3.deepcopy()
s7.fft_filter(
    transfer_function=w_low * w_high,
    function_domain="frequency",
    shift=True
)

_, ax = plt.subplots(figsize=(13, 6), ncols=2)
ax[0].imshow(s3.inav[0, 0].data, cmap="gray")
ax[0].set_title("Static + dynamic background removed")
ax[1].imshow(s7.inav[0, 0].data, cmap="gray")
_ = ax[1].set_title("After FFT filtering")
FFT filtering:
[########################################] | 100% Completed |  0.3s
_images/pattern_processing_30_1.png

Note that filtering with a spatial kernel in the frequency domain, after creating the kernel’s transfer function via FFT, and computing the inverse FFT (IFFT), is, in this case, the same as spatially correlating the kernel with the pattern.

Let’s demonstrate this by attempting to sharpen a pattern with a Laplacian kernel in both the spatial and frequency domains and comparing the results (this is a purely illustrative example, and perhaps not that practically useful):

[15]:
s8 = s3.deepcopy()
w_laplacian = np.array([[-1, -1, -1], [-1, 8, -1], [-1, -1, -1]])
s8.fft_filter(transfer_function=w_laplacian, function_domain="spatial")
FFT filtering:
[########################################] | 100% Completed |  0.1s
[16]:
from scipy.ndimage import correlate


p_filt = correlate(s3.inav[0, 0].data.astype(np.float32), weights=w_laplacian)
p_filt_resc = kp.pattern.rescale_intensity(p_filt, dtype_out=np.uint8)

_, ax = plt.subplots(figsize=(13, 6), ncols=2)
ax[0].imshow(s8.inav[0, 0].data, cmap="gray")
ax[0].set_title("Filter in frequency domain")
ax[1].imshow(p_filt_resc, cmap="gray")
ax[1].set_title("Filter in spatial domain")

np.sum(s8.inav[0, 0].data - p_filt_resc)  # Which is zero
[16]:
0
_images/pattern_processing_33_1.png

Note also that fft_filter() performs the filtering on the patterns with data type np.float32, and therefore have to rescale back to the pattern’s original data type if necessary.


Rescale intensity

Vendors usually write patterns to file with 8 (uint8) or 16 (uint16) bit integer depth, holding [0, 2\(^8\)] or [0, 2\(^{16}\)] gray levels, respectively. To avoid losing intensity information when processing, we often change data types to e.g. 32 bit floating point (float32). However, only changing the data type with change_dtype() does not rescale pattern intensities, leading to patterns not using the full available data type range:

[17]:
s9 = s3.deepcopy()

print(s9.data.dtype, s9.data.max())

s9.change_dtype(np.uint16)

print(s9.data.dtype, s9.data.max())

plt.figure(figsize=(6, 5))
plt.imshow(s9.inav[0, 0].data, vmax=1000, cmap="gray")
plt.title("16-bit pixels w/ 255 as max intensity", pad=15)
_ = plt.colorbar()
uint8 255
uint16 255
_images/pattern_processing_36_1.png

In these cases it is convenient to rescale intensities to a desired data type range, either keeping relative intensities between patterns in a scan or not. We can do this for all patterns in an EBSD signal with kikuchipy.signals.EBSD.rescale_intensity():

[18]:
s9.rescale_intensity(relative=True)

print(s9.data.dtype, s9.data.max())

plt.figure(figsize=(6, 5))
plt.imshow(s9.inav[0, 0].data, cmap="gray")
plt.title("16-bit pixels w/ 65535 as max. intensity", pad=15)
_ = plt.colorbar()
Rescaling the image intensities:
[########################################] | 100% Completed |  0.4s
uint16 65535
_images/pattern_processing_38_1.png

Or, we can do it for a single pattern (numpy.ndarray) with kikuchipy.pattern.rescale_intensity():

[19]:
p = s3.inav[0, 0].data
print(p.min(), p.max())
p2 = kp.pattern.rescale_intensity(p, dtype_out=np.uint16)
print(p2.min(), p2.max())
0 255
0 65535

We can also stretch the pattern contrast by removing intensities outside a range passed to in_range or at certain percentiles by passing percentages to percentiles:

[20]:
s10 = s3.deepcopy()
print(s10.data.min(), s10.data.max())
s10.rescale_intensity(out_range=(10, 245))
print(s10.data.min(), s10.data.max())
0 255
Rescaling the image intensities:
[########################################] | 100% Completed |  0.1s
10 245
[21]:
s10.rescale_intensity(percentiles=(0.5, 99.5))
print(s10.data.min(), s3.data.max())
Rescaling the image intensities:
[########################################] | 100% Completed |  0.1s
0 255
[22]:
fig, ax = plt.subplots(figsize=(13, 5), ncols=2)
im0 = ax[0].imshow(s3.inav[0, 0].data, cmap="gray")
ax[0].set_title("Static + dynamic background removed", pad=15)
fig.colorbar(im0, ax=ax[0])
im1 = ax[1].imshow(s10.inav[0, 0].data, cmap="gray")
ax[1].set_title("After clipping lowest/highest intensities", pad=15)
_ = fig.colorbar(im1, ax=ax[1])
_images/pattern_processing_44_0.png

This can reduce the influence of outliers with exceptionally high or low intensities, like hot or dead pixels.


Normalize intensity

It can be useful to normalize pattern intensities to a mean value of \(\mu = 0.0\) and a standard deviation of e.g. \(\sigma = 1.0\) when e.g. comparing patterns or calculating the image quality. Patterns can be normalized with normalize_intensity():

[23]:
s11 = s3.deepcopy()
np.mean(s11.data)
[23]:
113.94771604938272
[24]:
# s11.change_dtype(np.float32)  # Or pass `dtype_out` as below
s11.normalize_intensity(num_std=1, dtype_out=np.float32)
np.mean(s11.data)
Normalizing the image intensities:
[########################################] | 100% Completed |  0.9s
[24]:
-1.9191225e-08
[25]:
_, ax = plt.subplots(figsize=(13, 4), ncols=2)
ax[0].hist(s3.data.ravel(), bins=255)
ax[0].set_title("Static + dynamic background removed")
ax[1].hist(s11.data.ravel(), bins=255)
_ = ax[1].set_title("After intensity normalization")
_images/pattern_processing_49_0.png