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

In [ ]:

```
# 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:

In [ ]:

```
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():

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().
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).

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`

:

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.

The Gaussian blurred pattern removed during dynamic background correction can be obtained as an EBSD signal by calling 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")
```

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():

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 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()
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):

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.

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:

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,
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`

:

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.

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:

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():

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():

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.

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():

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")
```