#!/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 matching
# 
# Crystal orientations can be determined from experimental EBSD patterns by
# matching them to simulated patterns of known phases and orientations, see
# e.g. <cite data-cite="chen2015dictionary">Chen et al. (2015)</cite>,
# <cite data-cite="nolze2016pattern">Nolze et al. (2016)</cite>, 
# <cite data-cite="foden2019indexing">Foden et al. (2019)</cite>.
# 
# Here, we will demonstrate *dictionary indexing* (DI) using a small Ni EBSD data
# set and a dynamically simulated Ni master pattern from EMsoft, both of
# low resolution and found in the [kikuchipy.data](../reference.rst#data) module.
# The pattern dictionary is generated from a uniform grid of orientations with a
# fixed projection center (PC). The true orientation is likely to fall in between
# grid points, which means there is always a lower angular accuracy associated
# with DI. We can improve upon each solution by letting the orientation deviate
# from the grid points. We do this by maximizing the similarity between
# experimental and simulated patterns using numerical optimization algorithms from
# the [SciPy library](https://docs.scipy.org/doc/scipy/reference/optimize.html).
# This is here called *orientation refinement*. We could instead keep the
# orientations fixed and let the PC parameters deviate from their fixed values used
# in the dictionary, here called *projection center refinement*. Finally, we can
# also refine both at the same time, here called *orientation and projection center
# refinement*. The need for orientation or orientation and PC refinement is
# discussed by e.g.
# <cite data-cite="singh2017application">Singh et al. (2017)</cite>,
# <cite data-cite="winkelmann2020refined">Winkelmann et al. (2020)</cite>, and
# <cite data-cite="pang2020optimization">Pang et al. (2020)</cite>.
# 
# The term *pattern matching* is here used for the combined approach of DI
# followed by refinement.
# 
# Before we can generate a dictionary of
# simulated patterns, we need a master pattern containing all possible scattering
# vectors for a candidate phase. This can be simulated using EMsoft
# (<cite data-cite="callahan2013dynamical">Callahan and De Graef (2013)</cite> and
# <cite data-cite="jackson2014h5ebsd">Jackson et al. (2014)</cite>) and then read
# into kikuchipy.

# First, we import libraries and load the small experimental Nickel test data

# In[1]:


# Exchange inline for notebook or qt5 (from pyqt) for interactive plotting
get_ipython().run_line_magic('matplotlib', 'inline')

import tempfile
import matplotlib.pyplot as plt
import hyperspy.api as hs
import numpy as np
from orix import sampling, plot, io
from orix.vector import Vector3d
import kikuchipy as kp


plt.rcParams.update({"figure.facecolor": "w", "font.size": 15})

# Use kp.load("data.h5") to load your own data
s = kp.data.nickel_ebsd_large(allow_download=True)  # External download
s


# To obtain a good match, we must increase the signal-to-noise ratio. In this
# pattern matching analysis, the Kikuchi bands are considered the signal, and the
# angle-dependent backscatter intensity, along with unwanted detector effects,
# are considered to be noise. See the
# [pattern processing guide](pattern_processing.rst) for further details.

# In[2]:


s.remove_static_background()
s.remove_dynamic_background()


# <div class="alert alert-info">
# 
# Note
# 
# DI is computationally intensive and takes in general a long time to run due to
# all the pattern comparisons being done. To maintain an acceptable memory usage
# and be done within reasonable time, it is recommended to write processed
# patterns to an HDF5 file for quick reading during DI.
# 
# </div>

# In[3]:


#s.save("pattern_static_dynamic.h5")


# ## Dictionary indexing
# 
# ### Load a master pattern
# 
# Next, we load a dynamically simulated Nickel master pattern generated with
# EMsoft, in the northern hemisphere projection of the square Lambert projection
# for an accelerating voltage of 20 keV

# In[4]:


energy = 20
mp = kp.data.nickel_ebsd_master_pattern_small(projection="lambert", energy=energy)
mp


# In[5]:


mp.plot()


# The Nickel phase information, specifically the crystal symmetry, asymmetric atom
# positions, and crystal lattice, is conveniently stored in an
# [orix.crystal_map.Phase](https://orix.readthedocs.io/en/stable/reference.html#orix.crystal_map.Phase)

# In[6]:


ni = mp.phase
ni


# In[7]:


ni.structure  # Element, x, y, z, site occupation


# In[8]:


ni.structure.lattice  # nm and degrees


# ### Sample orientation space
# 
# If we don't know anything about the possible crystal (unit cell) orientations in
# our sample, the safest thing to do is to generate a dictionary of orientations
# uniformly distributed in a candidate phase's orientation space. To achieve this,
# we sample the Rodrigues Fundamental Zone of the proper point group *432* using
# cubochoric sampling <cite data-cite="singh2016orientation">Singh and De Graef
# (2016)</cite> available in
# [orix.sampling.get_sample_fundamental()](https://orix.readthedocs.io/en/stable/reference.html#orix.sampling.get_sample_fundamental). We can
# choose the average disorientation between orientations, the "resolution", of
# this sampling. Here, we will use a rather low resolution of 3$^{\circ}$.
# 
# <div class="alert alert-info">
# 
# Note
# 
# Cubochoric sampling became available in orix v0.7.
# 
# </div>

# In[9]:


rotations = sampling.get_sample_fundamental(
    method="cubochoric", resolution=3, point_group=ni.point_group
)
rotations


# This sampling resulted in 30 443 crystal orientations. See the [orix user guide
# on orientation sampling](https://orix.readthedocs.io/en/stable/uniform_sampling_of_orientation_space.html)
# for further details and options for orientation sampling.
# 
# <div class="alert alert-info">
# 
# Note
# 
# An average disorientation of 3$^{\circ}$ results in a course sampling of
# orientation space; a lower average disorientation should be used for
# experimental work.
# 
# </div>

# ### Define the detector-sample-geometry

# Now that we have our master pattern and crystal orientations, we need to
# describe the EBSD detector's position with respect to the sample (interaction
# volume). This ensures that projecting parts of the master pattern onto our
# detector yields dynamically simulated patterns resembling our experimental ones.
# See the [reference frames](reference_frames.rst) user guide and the
# [EBSDDetector](../reference.rst#kikuchipy.detectors.EBSDDetector)
# class for further details.

# In[10]:


signal_shape = s.axes_manager.signal_shape[::-1]
detector = kp.detectors.EBSDDetector(
    shape=signal_shape,
    pc=[0.421, 0.7794, 0.5049],
    sample_tilt=70,
    convention="tsl",
)
detector


# Let's double check the projection/pattern center (PC) position on the detector
# using
# [plot()](../reference.rst#kikuchipy.detectors.EBSDDetector.plot)

# In[11]:


detector.plot(coordinates="gnomonic", pattern=s.inav[0, 0].data)


# ### Generate dictionary
# 
# Now we're ready to generate our dictionary of simulated patterns by projecting
# parts of the master pattern onto our detector for all sampled orientations,
# using the
# [get_patterns()](../reference.rst#kikuchipy.signals.ebsdmasterpattern.get_patterns)
# method. The method assumes the crystal orientations are represented with respect
# to the EDAX TSL sample reference frame RD-TD-ND.
# 
# <div class="alert alert-info">
# 
# Note
# 
# It is in general adviced to not compute the dictionary immediately, but let the
# dictionary indexing method handle this, by passing `compute=False`. This will
# return a `LazyEBSD` signal, with the dictionary patterns as a
# [dask array](https://docs.dask.org/en/latest/array.html).
# 
# </div>

# In[12]:


sim = mp.get_patterns(
    rotations=rotations,
    detector=detector,
    energy=energy,
    dtype_out=np.float32,
    compute=True
)
sim


# Let's inspect the three first of the 30 443 simulated patterns

# In[13]:


#sim.plot()  # Plot the patterns with a navigator for easy inspection
fig, ax = plt.subplots(ncols=3, figsize=(18, 6))
for i in range(3):
    ax[i].imshow(sim.inav[i].data, cmap="gray")
    euler = np.rad2deg(sim.xmap[i].rotations.to_euler())[0]
    ax[i].set_title(
        f"($\phi_1, \Phi, \phi_2)$ = {np.array_str(euler, precision=1)}"
    )
    ax[i].axis("off")
fig.tight_layout()


# ### Perform indexing
# 
# *Signal masking was added in version: 0.5.*
# 
# The Kikuchi pattern signal is usually weak towards the corners of the detector,
# so we can pass a signal mask to only match the pixels where the mask values are
# ``False``, i.e. mask out the values that are ``True``. This convention is in
# line with how NumPy, Dask, scikit-image etc. defines a mask. We can pass
# whatever mask we want, as long as it is a boolean array of the detector shape.

# In[14]:


signal_mask = ~kp.filters.Window("circular", signal_shape).astype(bool)

p = s.inav[0, 0].data
fig, ax = plt.subplots(ncols=2, figsize=(10, 5))
ax[0].imshow(p * signal_mask, cmap="gray")
ax[0].set_title("Not used in matching")
ax[1].imshow(p * ~signal_mask, cmap="gray")
ax[1].set_title("Used in matching");


# Finally, let's use the
# [dictionary_indexing()](../reference.rst#kikuchipy.signals.EBSD.dictionary_indexing)
# method to match the simulated patterns to our experimental patterns, using
# the [zero-mean normalized cross correlation (NCC)](../reference.rst#kikuchipy.indexing.similarity_metrics.NormalizedCrossCorrelationMetric)
# coefficient $r$
# <cite data-cite="gonzalez2017digital">Gonzalez & Woods (2017)</cite>, which is
# the default similarity metric. Let's keep the 20 best matching orientations. A
# number of 4125 * 30443 comparisons is quite small, which we can do in memory all
# at once. However, in cases where the number of comparisons are too big for our
# memory to handle, we should iterate over the dictionary of simulated patterns
# by passing the number of patterns per iteration. To demonstrate this, we do at
# least 10 iterations here. The results are returned as a
# [orix.crystal_map.CrystalMap](https://orix.readthedocs.io/en/stable/reference.html#crystalmap).

# In[15]:


xmap = s.dictionary_indexing(
    sim,
    metric="ncc",
    keep_n=20,
    n_per_iteration=sim.axes_manager.navigation_size // 10,
    signal_mask=signal_mask
)
xmap


# <div class="alert alert-info">
# 
# Note
# 
# *Added in version: 0.5.*
#     
# Dictionary indexing of real world data sets takes a long time because the
# matching is computationally intensive. The
# [dictionary_indexing()](../reference.rst#kikuchipy.signals.EBSD.dictionary_indexing)
# method accepts parameters *n_per_iteration*, *rechunk* and *dtype* that lets us
# control this behaviour to a certain extent, so be sure to take a look at the
# method's docstring.
# 
# </div>

# The [normalized dot product](../reference.rst#kikuchipy.indexing.similarity_metrics.NormalizedDotProductMetric)
# can be used instead of the NCC by passing `metric="ndp"`. A custom metric can
# be used instead, by creating a class which inherits from the abstract class
# [SimilarityMetric](../reference.rst#kikuchipy.indexing.similarity_metrics.SimilarityMetric).

# The results can be exported to an HDF5 file re-readable by orix, or an
# .ang file readable by MTEX and some commercial packages

# In[16]:


temp_dir = tempfile.mkdtemp() + "/"
io.save(temp_dir + "ni.h5", xmap)
io.save(temp_dir + "ni.ang", xmap)


# ### Analyze indexing results
# 
# With the [orix library](https://orix.readthedocs.io) we can plot inverse pole
# figures (IPFs), color orientations to produce orientation maps (also called IPF
# maps), and more. If we want to inspect the results as grains, orientation
# density functions, or something else, we have to use other software, like MTEX.
# 
# Let's generate an IPF color key and plot it

# In[17]:


ckey = plot.IPFColorKeyTSL(ni.point_group)
print(ckey)
ckey.plot()


# With this color key we can color orientations according to which
# crystal directions the sample direction $Z$ points in in every
# map pixel

# In[18]:


xmap.scan_unit = "um"
xmap.plot(ckey.orientation2color(xmap.orientations), remove_padding=True)


# With a few more lines, we can plot maps for all three sample directions

# In[19]:


ori = xmap.orientations
directions = Vector3d(((1, 0, 0), (0, 1, 0), (0, 0, 1)))

fig, axes = plt.subplots(figsize=(15, 5), ncols=3)
for ax, v, title in zip(axes, directions, ("X", "Y", "Z")):
    ckey.direction = v
    rgb = ckey.orientation2color(ori).reshape(xmap.shape + (3,))
    ax.imshow(rgb)
    ax.axis("off")
    ax.set_title(f"IPF {title}")
fig.tight_layout()


# The sample is recrystallized Ni, so we expect a continuous color within
# grains, except for twinning grains. The orientation maps are mostly in
# line with our expectation. Some grains have a scatter of similar
# colors, which is most likely due to the discrete nature of our
# dictionary. To improve upon this result we can reduce the orientation
# sampling characteristic distance. This will increase our dictionary size
# and thus our indexing time. An alternative is to perform orientation
# refinement, which we will do below.
# 
# We can get further confirmation of our initial analysis by inspecting
# some indexing quality maps, like the best matching scores and so-called
# orientation similarity (OS) map, which compares the best matching
# orientations for each pattern to it's nearest neighbours. Let's get the
# NCC map in the correct shape from the CrystalMap’s scores property and
# the OS map with
# [orientation_similarity_map()](../reference.rst#kikuchipy.indexing.orientation_similarity_map)
# using the full list of best matches per point

# In[20]:


ncc_map = xmap.scores[:, 0].reshape(*xmap.shape)
os_map = kp.indexing.orientation_similarity_map(xmap)


# And plot the maps

# In[21]:


fig, ax = plt.subplots(ncols=2, figsize=(10, 3))
im0 = ax[0].imshow(ncc_map)
im1 = ax[1].imshow(os_map)
fig.colorbar(im0, ax=ax[0], label="NCC")
fig.colorbar(im1, ax=ax[1], label="Orientation similarity")
ax[0].axis("off")
ax[1].axis("off")
fig.tight_layout()


# From the NCC map we see that some grains match better than others. Due to
# the discrete nature of our dictionary, it is safe to assume that the best
# matching grains have patterns more similar to those in the dictionary than
# the others, that is, the different coefficients does not reflect anything
# physical in the microstructure. The OS map doesn't tell us much more than
# the NCC map in this case.
# 
# We can use the crystal map property `simulation_indices` to get the best
# matching simulated patterns from the dictionary of simulated patterns

# In[22]:


best_patterns = sim.data[xmap.simulation_indices[:, 0]].reshape(s.data.shape)
s_best = kp.signals.EBSD(best_patterns)
s_best


# The simplest way to visually compare the experimental and best matching
# simulated patterns are to
# [plot them in the same navigator](visualizing_patterns.ipynb#plot-multiple-signals).
# Let's create an RGB navigator signal from the IPF $Z$ orientation map with the
# convenience function
# [get_rgb_navigator()](../reference.rst#kikuchipy.draw.get_rgb_navigator). When
# using an interactive backend like `Qt5Agg`, we can then move the red square
# around to look at the patterns in each point.

# In[23]:


rgb_navigator = kp.draw.get_rgb_navigator(rgb)
hs.plot.plot_signals([s, s_best], navigator=rgb_navigator)


# Let's also plot the best matches for patterns from two grains

# In[24]:


grain1 = (0, 0)
grain2 = (30, 10)
fig, ax = plt.subplots(ncols=2, nrows=2, figsize=(10, 10))
ax[0, 0].imshow(s.inav[grain1].data, cmap="gray")
ax[0, 0].axis("off")
ax[0, 1].imshow(s_best.inav[grain1].data, cmap="gray")
ax[0, 1].axis("off")
ax[1, 0].imshow(s.inav[grain2].data, cmap="gray")
ax[1, 0].axis("off")
ax[1, 1].imshow(s_best.inav[grain2].data, cmap="gray")
ax[1, 1].axis("off")
fig.tight_layout(h_pad=0.5, w_pad=1)


# ## Refinement
# 
# Let's look at the effect of three refinement routes, all implemented as `EBSD`
# methods:
# 
# 1. Refine orientations while keeping the PCs fixed: [refine_orientation()](../reference.rst#kikuchipy.signals.EBSD.refine_orientation) 
# 2. Refine PCs while keeping the orientations fixed: [refine_projection_center()](../reference.rst#kikuchipy.signals.EBSD.refine_projection_center) 
# 3. Refine orientations and PCs at the same time: [refine_orientation_projection_center()](../reference.rst#kikuchipy.signals.EBSD.refine_orientation_projection_center) 
# 
# For each run we will compare the maps and histograms of NCC scores before and
# after refinement, and also the PC parameters if appropriate.
# 
# As stated at the top, we use the numerical optimization routines from the
# [SciPy library](https://docs.scipy.org/doc/scipy/reference/optimize.html). For
# every orientation and/or PC, we want to iteratively increase the similarity (NCC
# score) between our experimental pattern and a new simulated pattern projected
# onto our EBSD detector for every iteration until the score increase from one
# iteration to the next is below a certain threshold, by default 0.0001 in most
# cases. The orientation and/or PC is updated slightly for every iteration. We
# have access to both local and global optimization algorithms. Consult the
# kikuchipy docstring methods linked below and the SciPy documentation for all
# available options.
# 
# Note that while we here refine orientations obtained from DI, any orientation
# results, e.g. from Hough indexing, can be refined with this approach, as long as
# an appropriate master pattern and
# [EBSDDetector](../reference.rst#kikuchipy.detectors.EBSDDetector) is provided,
# and the orientation results are passed as a
# [CrystalMap](https://orix.readthedocs.io/en/stable/reference.html#orix.crystal_map.CrystalMap).
# 
# ### Refine orientations
# 
# We use
# [refine_orientation()](../reference.rst#kikuchipy.signals.EBSD.refine_orientation)
# to refine orientations while keeping the PCs fixed, using the default local
# Nelder-Mead simplex method

# In[25]:


xmap_refined = s.refine_orientation(
    xmap=xmap,
    detector=detector,
    master_pattern=mp,
    energy=energy,
    method="minimize",  # Default
    method_kwargs=dict(method="Nelder-Mead", options=dict(fatol=1e-4)),  # Default
    compute=True,  # Default
)


# Compare the NCC score maps. We use the same colorbar bounds for both maps

# In[26]:


ncc_after_orientation_refinement = xmap_refined.get_map_data("scores")

ncc_di_min = np.min(ncc_map)
ncc_di_max = np.max(ncc_map)
ncc_ori_ref_min = np.min(ncc_after_orientation_refinement)
ncc_ori_ref_max = np.max(ncc_after_orientation_refinement)

vmin = min([ncc_di_min, ncc_ori_ref_min])
vmax = max([ncc_di_max, ncc_ori_ref_max])


# In[27]:


ncc_after_label = "NCC after ori. ref."

fig, ax = plt.subplots(ncols=2, figsize=(10, 3))
im0 = ax[0].imshow(ncc_map, vmin=vmin, vmax=vmax)
im1 = ax[1].imshow(ncc_after_orientation_refinement, vmin=vmin, vmax=vmax)
fig.colorbar(im0, ax=ax[0], label="NCC from DI")
fig.colorbar(im1, ax=ax[1], label=ncc_after_label)
for a in ax:
    a.axis("off")
fig.tight_layout(w_pad=-1);


# Compare the histograms

# In[28]:


bins = np.linspace(vmin, vmax, 100)
fig, ax = plt.subplots(figsize=(9, 5))
_ = ax.hist(ncc_map.ravel(), bins, alpha=0.5, label="NCC from DI")
_ = ax.hist(
    ncc_after_orientation_refinement.ravel(),
    bins,
    alpha=0.5,
    label=ncc_after_label,
)
ax.set_xlabel("Normalized cross correlation (NCC) scores")
ax.set_ylabel("Frequency")
ax.legend()
fig.tight_layout();


# We see that DI found the best orientation (with a fixed PC) for most of the
# patterns, which the refinement was able to improve further. However, there
# are a few patterns with a very low NCC score (0.1-0.2) which refinement
# couldn't improve upon, which tells us that these were misindexed during DI.
# If there are Kikuchi bands in these patterns, a larger dictionary with a
# finer orientation sampling could improve indexing of them.
# 
# Let's also plot the (IPF) $X$ orientation maps before and after refinement, and
# also the disorientation angle (smallest misorientation angle found after
# accounting for symmetry) difference between the two maps

# In[29]:


ori_ref = xmap_refined.orientations

ckey.direction = Vector3d.xvector()
rgb_z = ckey.orientation2color(ori)
rgb_z_ref = ckey.orientation2color(ori_ref)

mori_angle = np.rad2deg(ori.angle_with(ori_ref).data)

fig, ax = plt.subplots(ncols=3, figsize=(14, 3))
ax[0].imshow(rgb_z.reshape(xmap.shape + (3,)))
ax[1].imshow(rgb_z_ref.reshape(xmap.shape + (3,)))
im2 = ax[2].imshow(mori_angle.reshape(xmap.shape))
fig.colorbar(im2, ax=ax[2], label=r"Disorientation angle $\omega$ [$^{\circ}$]")
for a in ax:
    a.axis("off")
fig.tight_layout(w_pad=-10)


# We see that refinement was able to remove the scatter of similar colors
# in the grains.

# ### Refine projection centers
# 
# We use
# [refine_projection_center()](../reference.rst#kikuchipy.signals.EBSD.refine_projection_center)
# to refine PCs while keeping the orientations fixed, using the local modified
# Powell method. This method is also known as Bound Optimization BY Quadratic
# Approximation (BOBYQA), and is used in EMsoft and discussed by
# <cite data-cite="singh2017application">Singh et al. (2017)</cite>. We will pass
# a `trust_region` of +/- 5% for the PC parameters to use as upper and lower
# bounds during refinement.
# 
# We can also pass `compute=False`, to do the
# refinement later. The results will then be a
# [dask.array.Array](https://docs.dask.org/en/latest/generated/dask.array.Array.html#dask.array.Array).
# We can pass this array to
# [kikuchipy.indexing.compute_refine_projection_center_results()](../reference.rst#kikuchipy.indexing.compute_refine_projection_center_results)
# and perform the refinement to retrieve the results

# In[30]:


result_array = s.refine_projection_center(
    xmap=xmap,
    detector=detector,
    master_pattern=mp,
    energy=energy,
    method="minimize",
    method_kwargs=dict(method="Powell", options=dict(ftol=1e-3)),
    trust_region=[0.05, 0.05, 0.05],
    compute=False,
)


# In[31]:


ncc_after_pc_refinement, detector_refined = kp.indexing.compute_refine_projection_center_results(
    results=result_array, detector=detector, xmap=xmap
)


# Note that `refine_orientation()` and `refine_orientation_projection_center()`
# also takes the `compute` parameter, and there are similar functions to compute the
# refinement results:
# [kikuchipy.indexing.compute_refine_orientation_results()](../reference.rst#kikuchipy.indexing.compute_refine_orientation_results)
# and
# [kikuchipy.indexing.compute_refine_orientation_projection_center_results()](../reference.rst#kikuchipy.indexing.compute_refine_orientation_projection_center_results).
# 
# Let's plot the refined scores and PCs

# In[32]:


ncc_pc_ref_min = np.min(ncc_after_pc_refinement)
ncc_pc_ref_max = np.max(ncc_after_pc_refinement)

vmin2 = min([ncc_di_min, ncc_pc_ref_min])
vmax2 = max([ncc_di_max, ncc_pc_ref_max])


# In[33]:


ncc_after_pc_label = "NCC after PC refinement"

fig, ax = plt.subplots(ncols=2, figsize=(10, 3))
im0 = ax[0].imshow(ncc_map, vmin=vmin2, vmax=vmax2)
im1 = ax[1].imshow(ncc_after_pc_refinement, vmin=vmin2, vmax=vmax2)
ax[0].axis("off")
ax[1].axis("off")
fig.colorbar(im0, ax=ax[0], label="NCC from DI")
fig.colorbar(im1, ax=ax[1], label=ncc_after_pc_label)
fig.tight_layout();


# Compare the histograms

# In[34]:


bins = np.linspace(vmin2, vmax2, 100)
fig, ax = plt.subplots(figsize=(9, 5))
_ = ax.hist(ncc_map.ravel(), bins, alpha=0.5, label="NCC from DI")
_ = ax.hist(
    ncc_after_pc_refinement.ravel(),
    bins,
    alpha=0.5,
    label=ncc_after_pc_label,
)
ax.set_xlabel("Normalized cross correlation (NCC) scores")
ax.set_ylabel("Frequency")
ax.legend()
fig.tight_layout();


# In[35]:


print(
    f"PC used in DI:\t\t{detector.pc_average}\n"
    f"PC after PC refinement:\t{detector_refined.pc_average}"
)

fig, ax = plt.subplots(ncols=3, figsize=(15, 3))
im0 = ax[0].imshow(detector_refined.pcx)
im1 = ax[1].imshow(detector_refined.pcy)
im2 = ax[2].imshow(detector_refined.pcz)
fig.colorbar(im1, ax=ax[1], label="Projection center y")
fig.colorbar(im0, ax=ax[0], label="Projection center x")
fig.colorbar(im2, ax=ax[2], label="Projection center z")
for a in ax:
    a.axis("off")
fig.tight_layout();


# ### Refine orientations and projection centers
# 
# We use
# [refine_orientation_projection_center()](../reference.rst#kikuchipy.signals.EBSD.refine_orientation_projection_center)
# to refine orientations and PCs at the same time. For the purpose of saving time
# and computation resources in this demonstration, we'll do this for the upper left
# quarter of the data set only

# In[36]:


nav_shape = s.axes_manager.navigation_shape[::-1]
y1, x1 = np.array(nav_shape) // 2
slices2d = (slice(0, y1), slice(0, x1))


# In[37]:


s2 = s.inav[slices2d[::-1]]  # HyperSpy flips the usual NumPy (row, column) to (column, row)
s2


# In[38]:


xmap2 = xmap[slices2d]
xmap2


# In[39]:


xmap_refined2, detector_refined2 = s2.refine_orientation_projection_center(
    xmap=xmap2,
    detector=detector,
    master_pattern=mp,
    energy=energy,
    method="minimize",
    method_kwargs=dict(options=dict(fatol=1e-3)),
    compute=True,
)


# In[40]:


xmap_refined2


# Compare the NCC score maps. We use the same colorbar bounds for both maps

# In[41]:


ncc_after_orientation_pc_refinement = xmap_refined2.get_map_data("scores")
ncc_map2 = ncc_map[slices2d]

ncc2_di_min = np.min(ncc_map2)
ncc2_di_max = np.max(ncc_map2)
ncc_ori_pc_ref_min = np.min(ncc_after_orientation_pc_refinement)
ncc_ori_pc_ref_max = np.max(ncc_after_orientation_pc_refinement)

vmin3 = min([ncc2_di_min, ncc_ori_pc_ref_min])
vmax3 = max([ncc2_di_max, ncc_ori_pc_ref_max])


# In[42]:


ncc_after_ori_pc_label = "NCC after ori./PC ref."

fig, ax = plt.subplots(ncols=2, figsize=(9, 3))
im0 = ax[0].imshow(ncc_map2, vmin=vmin, vmax=vmax)
im1 = ax[1].imshow(ncc_after_orientation_pc_refinement, vmin=vmin3, vmax=vmax3)
fig.colorbar(im0, ax=ax[0], label="NCC from DI")
fig.colorbar(im1, ax=ax[1], label=ncc_after_ori_pc_label)
ax[0].axis("off")
ax[1].axis("off")
fig.tight_layout();


# Compare the histograms

# In[43]:


bins = np.linspace(vmin3, vmax3, 100)
fig, ax = plt.subplots(figsize=(9, 5))
_ = ax.hist(ncc_map2.ravel(), bins, alpha=0.5, label="NCC from DI")
_ = ax.hist(
    ncc_after_orientation_pc_refinement.ravel(),
    bins,
    alpha=0.5,
    label=ncc_after_ori_pc_label,
)
ax.set_xlabel("Normalized cross correlation (NCC) scores")
ax.set_ylabel("Frequency")
ax.legend()
fig.tight_layout();


# Let's also inspect the refined PC parameters

# In[44]:


print(
    f"PC used in DI:\t\t{detector.pc_average}\n"
    f"PC after PC refinement:\t{detector_refined2.pc_average}"
)

pc_to_plot = (detector_refined2.pcx, detector_refined2.pcy, detector_refined2.pcz)
fig, ax = plt.subplots(ncols=3, figsize=(15, 3))
for a, to_plot, label in zip(ax, pc_to_plot, ("x", "y", "z")):
    im = a.imshow(to_plot)
    fig.colorbar(im, ax=a, label=f"Projection center {label}")
    a.axis("off")
fig.tight_layout();


# In[45]:


# Remove files written to disk in this user guide
import os
os.remove(temp_dir + "ni.h5")
os.remove(temp_dir + "ni.ang")
os.rmdir(temp_dir)