#!/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.
# # Feature maps
#
# This section details methods for extracting information from pattern
# intensities, called feature maps (for lack of a better description).
#
# Let's import the necessary libraries and a 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
# Use kp.load("data.h5") to load your own data
s = kp.data.nickel_ebsd_large(allow_download=True) # External download
s
# ## Image quality
#
# The image quality metric $Q$ presented by Krieger Lassen
# Lassen (1994)
# can be calculated for an [EBSD](../reference.rst#kikuchipy.signals.EBSD) signal
# with
# [get_image_quality()](../reference.rst#kikuchipy.signals.EBSD.get_image_quality),
# or, for a single pattern (`numpy.ndarray`), with
# [get_image_quality()](../reference.rst#kikuchipy.pattern.get_image_quality).
# Following the notation in
# Marquardt et al. (2017), it
# is given by
#
# $$
# \begin{align}
# Q &= 1 - \frac{J}{J_{\mathrm{res}}w_{\mathrm{tot}}},\\
# J &= \sum_{h = -N/2}^{N/2} \sum_{k = -N/2}^{N/2} w(h, k)
# \left|\mathbf{q}\right|^2,\\
# J_{\mathrm{res}} &= \frac{1}{N^2} \sum_{h = -N/2}^{N/2}
# \sum_{k = -N/2}^{N/2} \left|\mathbf{q}\right|^2,\\
# w_{\mathrm{tot}} &= \sum_{h = -N/2}^{N/2} \sum_{k = -N/2}^{N/2} w(h, k).
# \end{align}
# $$
#
# The function $w(h, k)$ is the Fast Fourier Transform (FFT) power spectrum of the
# EBSD pattern, and the vectors $\mathbf{q}$ are the frequency vectors with
# components $(h, k)$. The sharper the Kikuchi bands, the greater the high
# frequency content of the power spectrum, and thus the closer $Q$ will be to
# unity.
#
# Since we want to use the image quality metric to determine how pronounced the
# Kikuchi bands in our patterns are, we first remove the static and dynamic
# background:
# In[ ]:
s.remove_static_background()
s.remove_dynamic_background()
# To visualize parts of the computation, we compute the power spectrum of a
# pattern in the Nickel EBSD data set and the frequency vectors, shift the
# zero-frequency components to the centre, and plot them:
# In[ ]:
p = s.inav[20, 11].data
p_fft = kp.pattern.fft(p, shift=True)
q = kp.pattern.fft_frequency_vectors(shape=p.shape)
fig, ax = plt.subplots(figsize=(22, 6), ncols=3)
title_kwargs = dict(fontsize=22, pad=15)
fig.subplots_adjust(wspace=0.05)
im0 = ax[0].imshow(p, cmap="gray")
ax[0].set_title("Pattern", **title_kwargs)
fig.colorbar(im0, ax=ax[0])
im1 = ax[1].imshow(np.log(kp.pattern.fft_spectrum(p_fft)), cmap="gray")
ax[1].set_title("Log of shifted power spectrum of FFT", **title_kwargs)
fig.colorbar(im1, ax=ax[1])
im2 = ax[2].imshow(np.fft.fftshift(q), cmap="gray")
ax[2].set_title(r"Shifted frequency vectors $q$", **title_kwargs)
_ = fig.colorbar(im2, ax=ax[2])
# If we don't want the EBSD patterns to be
# [zero-mean normalized](pattern_processing.ipynb#normalize-intensity) before
# computing $Q$, we must pass `get_image_quality(normalize=False)`.
#
# Let's compute the image quality $Q$ and plot it for the entire data set:
# In[ ]:
iq = s.get_image_quality(normalize=True) # Default
# In[ ]:
plt.figure(figsize=(10, 6))
plt.imshow(iq, cmap="gray")
plt.colorbar(label=r"Image quality, $Q$", pad=0.01)
_ = plt.axis("off")
# If we want to use this map to navigate around in when plotting patterns, we can
# easily do that as explained in the
# [visualizing patterns](visualizing_patterns.rst) guide.
# ## Average dot product
#
# The average dot product, or normalized cross-correlation when centering each
# pattern's intensity about zero and normalizing the intensities to a standard
# deviation $\sigma$ of 1 (which is the default behaviour), between each pattern
# and their four nearest neighbours, can be obtained for an
# [EBSD](../reference.rst#kikuchipy.signals.EBSD) signal with
# [get_average_neighbour_dot_product_map()](../reference.rst#kikuchipy.signals.EBSD.get_average_neighbour_dot_product_map)
# In[ ]:
adp = s.get_average_neighbour_dot_product_map()
# In[ ]:
plt.figure(figsize=(10, 6))
plt.imshow(adp, cmap="gray")
plt.colorbar(label="Average dot product", pad=0.01)
_ = plt.axis("off")
# The map displays how similar each pattern is to its neighbours. Grain
# boundaries, and some scratches on the sample, can be clearly seen as pixels with
# a lower value, signifying that they are more dissimilar to their neighbouring
# pixels, as the ones within grains where the neighbour pixel similarity is high.
#
# The map above was created by averaging the dot product matrix per map point,
# created by calculating the dot product between each pattern and their four
# nearest neighbours, which can be seen in the black spots (uneven sample surface)
# in the left grains
# In[ ]:
w1 = kp.filters.Window()
w1.plot()
# We could instead average with e.g. the eight nearest neighbours
# In[ ]:
w2 = kp.filters.Window(window="rectangular", shape=(3, 3))
w2.plot()
# In[ ]:
adp2 = s.get_average_neighbour_dot_product_map(window=w2)
plt.figure(figsize=(10, 6))
plt.imshow(adp2, cmap="gray")
plt.colorbar(label="Average dot product", pad=0.01)
_ = plt.axis("off")
# Note that the window coefficients must be integers.
#
# We can also control whether pattern intensities should be centered about zero
# and/or whether they should be normalized prior to calculating the dot products
# by passing `zero_mean=False` and/or `normalize=False`. These are `True` by
# default. The data type of the output map, 32-bit floating point by default,
# can be set by passing e.g. `dtype_out=np.float64`.
#
# We can obtain the dot product matrices per map point, that is the matrices
# before they are averaged, with
# [get_neighbour_dot_product_matrices()](../reference.rst#kikuchipy.signals.EBSD.get_neighbour_dot_product_matrices).
# Let's see similar a pattern on a grain boundary in map point (x, y) = (50, 19)
# is to all its nearest neighbour in a (5, 5) window centered on that point
# In[ ]:
w3 = kp.filters.Window("rectangular", shape=(5, 5))
dp_matrices = s.get_neighbour_dot_product_matrices(window=w3)
# In[ ]:
x, y = (50, 19)
s_dp_matrices = hs.signals.Signal2D(dp_matrices)
s_dp_matrices.inav[x, y].plot()
# We can see that the pattern is more similar to the patterns up to the right,
# while it is quite dissimilar to the patterns to the lower left. Let's visualize
# this more clearly, as is done e.g. in Fig. 1 by
# Brewick et al. (2019)
# In[ ]:
y_n, x_n = w3.n_neighbours
s2 = s.inav[x - x_n : x + x_n + 1, y - y_n : y + y_n + 1].deepcopy()
s2.rescale_intensity(percentiles=(0.5, 99.5)) # Stretch the contrast a bit
s3 = s2 * s_dp_matrices.inav[x, y].T # Signals must have same navigation shape
# In[ ]:
_ = hs.plot.plot_images(
images=s3,
per_row=5,
label=None,
suptitle=None,
axes_decor=None,
colorbar=None,
vmin=int(s3.data.min()),
vmax=int(s3.data.max()),
padding=dict(wspace=0, hspace=-0.05),
fig=plt.figure(figsize=(10, 10)),
)
# Finally, we can pass this dot product matrix directly to
# [get_average_neighbour_dot_product_map()](../reference.rst#kikuchipy.signals.EBSD.get_average_neighbour_dot_product_map)
# via the `dp_matrices` parameter to obtain the average dot product map from these
# matrices
# In[ ]:
adp3 = s.get_average_neighbour_dot_product_map(dp_matrices=dp_matrices)
# Let's plot this and highlight the location of the pattern on the grain boundary
# above with a red circle
# In[ ]:
plt.figure(figsize=(10, 6))
plt.imshow(adp3, cmap="gray")
plt.colorbar(label="Average dot product", pad=0.01)
plt.scatter(x=x, y=y, marker="o", c="r", s=50)
_ = plt.axis("off")