#!/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 qt5 for interactive plotting from the pyqt package
get_ipython().run_line_magic('matplotlib', 'inline')

import hyperspy.api as hs
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


# 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].
# 
# 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)
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")


# ---
# ## 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
# <cite data-cite="gonzalez2017digital">Gonzalez and Woods (2017)</cite>
# 
# \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.Window) object, can be plotted:

# In[ ]:


w = kp.filters.Window(window="gaussian", shape=(3, 3), std=1)
_, _, _ = w.plot(cmap="inferno")


# 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.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.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[ ]:


fig, im, cbar = w.plot()


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

# <div class="alert alert-info">
# 
# 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
# <cite data-cite="wright2015introduction">Wright et al. (2015)</cite> for a
# discussion of this concern.
# 
# </div>

# ---
# ## Adaptive histogram equalization
# 
# Enhancing the pattern contrast with adaptive histogram equalization has been
# found useful when comparing patterns for dictionary indexing
# <cite data-cite="marquardt2017quantitative">Marquardt et al. (2017)</cite>.
# 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
# <cite data-cite="jackson2019dictionary">Jackson et al. (2019)</cite>, 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.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")