#!/usr/bin/env python
# coding: utf-8
# This notebook is part of the `kikuchipy` documentation https://kikuchipy.org.
# Links to the documentation won't work from the notebook.
# # 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](virtual_backscatter_electron_imaging.rst),
# 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](../reference.rst#kikuchipy.signals.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](../reference.rst#pattern) module.
#
# Let's import the necessary libraries and read the Nickel EBSD test data set
# In[ ]:
# Exchange inline for notebook or qt5 (from pyqt) for interactive plotting
get_ipython().run_line_magic('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()](http://hyperspy.org/hyperspy-doc/current/api/hyperspy.signal.html#hyperspy.signal.BaseSignal.deepcopy)
# of the original signal. As an example here, we create a new EBSD signal from a
# small part of the original signal:
# In[ ]:
s2 = s.deepcopy()
np.may_share_memory(s.data, s2.data)
# ---
#
# ## 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()](../reference.rst#kikuchipy.signals.EBSD.remove_static_background):
# In[ ]:
s2.remove_static_background(operation="subtract", relative=True)
# In[ ]:
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 as part of the
# signal `metadata`, which can be loaded via
# [set_experimental_parameters()](../reference.rst#kikuchipy.signals.EBSD.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].
#
#
#
# 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](../reference.rst#kikuchipy.signals.EBSD.xmap) and
# [EBSD.detector](../reference.rst#kikuchipy.signals.EBSD.detector) properties,
# which keeps track of the
# [CrystalMap](https://orix.readthedocs.io/en/stable/reference.html#orix.crystal_map.CrystalMap)
# and [EBSDDetector](../reference.rst#kikuchipy.detectors.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`, 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()](../reference.rst#kikuchipy.signals.EBSD.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](#Filtering-in-the-frequency-domain). The
# individual Gaussian blurred dynamic backgrounds are then subtracted or divided
# from the respective patterns, set by `operation`:
# In[ ]:
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.
# ---
# ## 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()](../reference.rst#kikuchipy.signals.EBSD.get_dynamic_background):
# In[ ]:
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")
# ---
# ## 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()](../reference.rst#kikuchipy.signals.EBSD.average_neighbour_patterns):
# In[ ]:
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](../reference.rst#kikuchipy.filters.Window) object, can be plotted:
# In[ ]:
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()](../reference.rst#kikuchipy.signals.EBSD.average_neighbour_patterns)
# or [Window](../reference.rst#kikuchipy.filters.Window) as a `numpy.ndarray`
# or a `dask.array.Array`. Additionally, any window in
# [scipy.signal.windows.get_window()](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.get_window.html)
# 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()](../reference.rst#kikuchipy.filters.Window.make_circular),
# although we could create a circular window directly by calling
# `window="circular"` upon window initialization):
# In[ ]:
w = kp.filters.Window(window="rectangular", shape=(5, 4))
w
# In[ ]:
w.make_circular()
w
# In[ ]:
s5 = s3.deepcopy()
s5.average_neighbour_patterns(w)
# In[ ]:
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.
#
#
# ---
# ## Adaptive histogram equalization
#
# 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()](../reference.rst#kikuchipy.signals.EBSD.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:
# In[ ]:
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](https://github.com/EMsoft-org/EMsoft/wiki/DItutorial#52-determination-of-pattern-pre-processing-parameters),
# 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()](../reference.rst#kikuchipy.signals.EBSD.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](../reference.rst#kikuchipy.filters.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`:
# In[ ]:
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:
# In[ ]:
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):
# In[ ]:
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")
# In[ ]:
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.
# ---
# ## 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()](http://hyperspy.org/hyperspy-doc/current/api/hyperspy.signal.html#hyperspy.signal.BaseSignal.change_dtype) does not
# rescale pattern intensities, leading to patterns not using the full available
# data type range:
# In[ ]:
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()](../reference.rst#kikuchipy.signals.EBSD.rescale_intensity):
# In[ ]:
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()](../reference.rst#kikuchipy.pattern.rescale_intensity):
# In[ ]:
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`:
# In[ ]:
s10 = s3.deepcopy()
print(s10.data.min(), s10.data.max())
s10.rescale_intensity(out_range=(10, 245))
print(s10.data.min(), s10.data.max())
# In[ ]:
s10.rescale_intensity(percentiles=(0.5, 99.5))
print(s10.data.min(), s3.data.max())
# In[ ]:
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.
# ---
# ## 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](feature_maps.ipynb#image-quality). Patterns can be
# normalized with
# [normalize_intensity()](../reference.rst#kikuchipy.signals.EBSD.normalize_intensity):
# In[ ]:
s11 = s3.deepcopy()
np.mean(s11.data)
# In[ ]:
# 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)
# In[ ]:
_, 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")