This notebook is part of the kikuchipy
documentation https://kikuchipy.org.
Links to the documentation won't work from the notebook.
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
# Exchange inline for notebook or qt5 (from pyqt) for interactive plotting
%matplotlib inline
import hyperspy.api as hs
import matplotlib.pyplot as plt
import numpy as np
import kikuchipy as kp
s = kp.data.nickel_ebsd_small() # Use kp.load("data.h5") to load your own data
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:
s2 = s.deepcopy()
np.may_share_memory(s.data, s2.data)
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():
s2.remove_static_background()
fig, 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[0].axis("off")
ax[1].imshow(s2.inav[0, 0].data, cmap="gray")
ax[1].set_title("Static background removed")
ax[1].axis("off")
fig.tight_layout(w_pad=0 - 8)
Here, the static background pattern is assumed to be stored in the
EBSD.static_background
property, which can be set and changed as desired. The static background pattern
can also be passed to the static_bg
parameter.
Since version 0.6, relative intensities between patterns cannot be maintained
via the relative
parameter. This was because the initial implementation filled
the RAM when removing the background from a large dataset, and a way around this
was not found.
Warning
The Acquisition_instrument.SEM.Detector.EBSD
and Sample.Phases
metadata
nodes are deprecated and will be removed in v0.6.
There are three main reasons for this change: the first is that only the static
background array stored in the
Acquisition_instrument.SEM.Detector.EBSD.static_background
node is used
internally, and so the remaining metadata is unnecessary. The background pattern
will be stored in its own EBSD.static_background
property instead. The second
is that keeping track of the unnecessary metadata makes writing and maintaining
input/ouput plugins challenging. The third is that the
EBSD.xmap and
EBSD.detector properties,
which keeps track of the
CrystalMap
and EBSDDetector for a
signal, respectively, should be used instead of the more "static" metadata.
The static background pattern intensities can be rescaled to each individual
pattern's intensity range before removal by passing scale_bg=True
.
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
:
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
)
# _ means we don't want the output
_, 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")
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.
The Gaussian blurred pattern removed during dynamic background correction can be obtained as an EBSD signal by calling get_dynamic_background():
bg = s.get_dynamic_background(filter_domain="frequency", std=8, truncate=4)
_, 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")
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():
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")
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 Gonzalez and Woods (2017)
\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:
w = kp.filters.Window(window="gaussian", shape=(3, 3), std=1)
w.plot()
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):
w = kp.filters.Window(window="rectangular", shape=(5, 4))
w
w.make_circular()
w
s5 = s3.deepcopy()
s5.average_neighbour_patterns(w)
w.plot()
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 Wright et al. (2015) for a discussion of this concern.
Enhancing the pattern contrast with adaptive histogram equalization has been
found useful when comparing patterns for dictionary indexing
Marquardt et al. (2017).
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:
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")
The kernel_size
parameter determines the size of the contextual regions. See
e.g. Fig. 5 in
Jackson et al. (2019), also
available via
EMsoft's GitHub repository wiki,
for the effect of varying kernel_size
.
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
:
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])
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:
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")
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):
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")
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
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.
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:
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()
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():
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()
Or, we can do it for a single pattern (numpy.ndarray
) with
kikuchipy.pattern.rescale_intensity():
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())
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
:
s10 = s3.deepcopy()
print(s10.data.min(), s10.data.max())
s10.rescale_intensity(out_range=(10, 245))
print(s10.data.min(), s10.data.max())
s10.rescale_intensity(percentiles=(0.5, 99.5))
print(s10.data.min(), s3.data.max())
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])
This can reduce the influence of outliers with exceptionally high or low intensities, like hot or dead pixels.
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():
s11 = s3.deepcopy()
np.mean(s11.data)
# 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)
_, 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")