{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "This notebook is part of the `kikuchipy` documentation https://kikuchipy.org.\n",
    "Links to the documentation won't work from the notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Geometrical EBSD simulations\n",
    "\n",
    "This section details how to inspect and visualize the results from pattern\n",
    "matching or Hough indexing **of cubic crystals** by plotting Kikuchi bands and\n",
    "zone axes onto an EBSD signal. We consider this a geometrical EBSD simulation,\n",
    "since it's only the Kikuchi band centres and zone axis positions that will be\n",
    "computed. These simulations are based on the work by Aimo Winkelmann in the\n",
    "supplementary material to\n",
    "<cite data-cite=\"britton2016tutorial\">Britton et al. (2016)</cite>.\n",
    "\n",
    "Let's import the necessary libraries and a small (3, 3) Nickel EBSD test data\n",
    "set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exchange inline for notebook or qt5 (from pyqt) for interactive plotting\n",
    "%matplotlib inline\n",
    "\n",
    "import tempfile\n",
    "from diffsims.crystallography import ReciprocalLatticePoint\n",
    "import hyperspy.api as hs\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from orix.crystal_map import Phase\n",
    "from orix.quaternion import Rotation\n",
    "import kikuchipy as kp\n",
    "\n",
    "\n",
    "plt.rcParams.update({\n",
    "    \"figure.figsize\": (10, 10), \"font.size\": 20, \"lines.markersize\": 10\n",
    "})\n",
    "\n",
    "s = kp.data.nickel_ebsd_small()  # Use kp.load(\"data.h5\") to load your own data\n",
    "s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's enhance the Kikuchi bands by removing the static and dynamic backgrounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s.remove_static_background()\n",
    "s.remove_dynamic_background()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = hs.plot.plot_images(\n",
    "    s, axes_decor=None, label=None, colorbar=False, tight_layout=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To project Kikuchi bands and zone axis positions onto our detector, we need\n",
    "1. a description of the crystal phase, in the geometrical case only the space\n",
    "   group\n",
    "2. the set of Kikuchi bands to consider, e.g. the {111}, {200}, {220}, and {311}\n",
    "   crystal plane families\n",
    "3. the crystal orientations with respect to the reference frame\n",
    "4. the position of the detector with respect to the sample reference frame,\n",
    "   in the form of a sample-detector model which includes the sample and detector\n",
    "   tilt and the projection center (shortes distance from the source point on the\n",
    "   sample to the detector), given here as (PC$_x$, PC$_y$, PC$_z$)\n",
    "\n",
    "We'll store the crystal phase information in an\n",
    "[orix.crystal_map.Phase](https://orix.readthedocs.io/en/stable/reference.html#orix.crystal_map.Phase)\n",
    "instance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phase = Phase(name=\"ni\", space_group=225)\n",
    "phase"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll set up the relevant Kikuchi bands (and the zone axes from these) in a\n",
    "[diffsims.crystallography.ReciprocalLatticePoint](https://diffsims.readthedocs.io/en/latest/reference.html#diffsims.crystallography.ReciprocalLatticePoint)\n",
    "instance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rlp = ReciprocalLatticePoint(\n",
    "    phase=phase, hkl=[[1, 1, 1], [2, 0, 0], [2, 2, 0], [3, 1, 1]]\n",
    ")\n",
    "rlp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot where the $(hkl)$ plane normals impinge on the unit sphere\n",
    "by showing their stereographic projection on the upper and lower\n",
    "hemispheres"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scatter_kwargs = dict(grid=True, hemisphere=\"both\", axes_labels=[\"x\", \"y\", None])\n",
    "colors = [f\"C{i}\" for i in range(rlp.size)]\n",
    "rlp.hkl.scatter(\n",
    "    c=colors,\n",
    "    vector_labels=[str(hi).replace(\"[\", \"(\").replace(\"]\", \")\") for hi in rlp.hkl.data],\n",
    "    **scatter_kwargs\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can get a new instance with the symmetrically equivalent planes $\\{hkl\\}$\n",
    "in each plane family using\n",
    "[ReciprocalLatticePoint.symmetrise()](https://diffsims.readthedocs.io/en/latest/reference.html#diffsims.crystallography.ReciprocalLatticePoint.symmetrise)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rlp2, mult = rlp.symmetrise(return_multiplicity=True)\n",
    "\n",
    "# Plot {111} plane family\n",
    "fig = rlp2.hkl[:mult[0]].scatter(c=colors[0], return_figure=True, **scatter_kwargs)\n",
    "\n",
    "# Plot remaining plane families\n",
    "i = 0\n",
    "for j in range(mult.size):\n",
    "    rlp2.hkl[i:i + mult[j]].scatter(c=colors[j], figure=fig)\n",
    "    i += mult[j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know from [pattern matching](pattern_matching.ipynb) of these nine patterns\n",
    "to dynamically simulated patterns of orientations uniformly\n",
    "distributed in the orientation space of the proper point group $432$, that they\n",
    "come from two grains with orientations of about $(\\phi_1, \\Phi, \\phi_2) =\n",
    "(80^{\\circ}, 34^{\\circ}, -90^{\\circ})$ and $(\\phi_1, \\Phi, \\phi_2) =\n",
    "(115^{\\circ}, 27^{\\circ}, -95^{\\circ})$. We store these orientations in an\n",
    "[orix.quaternion.Rotation](https://orix.readthedocs.io/en/stable/reference.html#orix.quaternion.Rotation)\n",
    "instance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grain1 = (80, 34, -90)\n",
    "grain2 = (115, 27, -95)\n",
    "rot = Rotation.from_euler(\n",
    "    np.radians([\n",
    "        [grain1, grain2, grain2],\n",
    "        [grain1, grain2, grain2],\n",
    "        [grain1, grain2, grain2]\n",
    "    ]),\n",
    ")\n",
    "rot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We describe the sample-detector model in an\n",
    "[kikuchipy.detectors.EBSDDetector](../reference.rst#kikuchipy.detectors.EBSDDetector)\n",
    "instance. From Hough indexing we know the projection center to be, in the EDAX\n",
    "TSL convention (see the [reference frame](reference_frames.rst) guide for the\n",
    "various conventions and more details on the use of the sample-detector model),\n",
    "$(x^{*}, y^{*}, z^{*}) = (0.421, 0.7794, 0.5049)$. The sample was tilted\n",
    "$70^{\\circ}$ about the microscope X direction towards the detector, and the\n",
    "detector normal was orthogonal to the optical axis (beam direction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "detector = kp.detectors.EBSDDetector(\n",
    "    shape=s.axes_manager.signal_shape[::-1],\n",
    "    sample_tilt=70,\n",
    "    pc=[0.421, 0.7794, 0.5049],\n",
    "    convention=\"tsl\"\n",
    ")\n",
    "detector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the projection center gets converted internally to the Bruker\n",
    "convention.\n",
    "\n",
    "Let's create an\n",
    "[EBSDSimulationGenerator](../reference.rst#kikuchipy.generators.EBSDSimulationGenerator)\n",
    "instance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_gen = kp.generators.EBSDSimulationGenerator(\n",
    "    detector=detector,\n",
    "    phase=phase,\n",
    "    rotations=rot\n",
    ")\n",
    "sim_gen"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're ready to simulate geometrical EBSD patterns from the generator and the\n",
    "sets of crystal plane families"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = sim_gen.geometrical_simulation(reciprocal_lattice_point=rlp2)\n",
    "sim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that 27 of the 50 Kikuchi bands we're visible on the detector in the nine\n",
    "patterns, stored in an instance of the\n",
    "[kikuchipy.simulations.features.KikuchiBand](../reference.rst#kikuchipy.simulations.features.KikuchiBand)\n",
    "class, which is a class inheriting from the\n",
    "[ReciprocalLatticePoint](https://diffsims.readthedocs.io/en/latest/reference.html#diffsims.crystallography.ReciprocalLatticePoint)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.zone_axes.size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also see that there are 91 zone axes resulting from the 27 Kikuchi bands,\n",
    "stored in an instance of the\n",
    "[kikuchipy.simulations.features.ZoneAxis](../reference.rst#kikuchipy.simulations.features.ZoneAxis)\n",
    "class, also inheriting from `ReciprocalLatticePoint`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now add these simulations as markers to be displayed on top of our\n",
    "experimental EBSD signal when plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "markers = sim.as_markers(pc=False)\n",
    "s.add_marker(marker=markers, plot_marker=False, permanent=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The markers update with the pattern when navigating, thus helping us determine\n",
    "whether an indexing was successful, and in labeling the bands and zone axes in\n",
    "the pattern"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s.plot(navigator=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Whether to plot only bands, zone axes, zone axes labels, projection center, or\n",
    "all of them, can be set in the\n",
    "[GeometricalEBSDSimulation.as_markers()](../reference.rst#kikuchipy.simulations.GeometricalEBSDSimulation.as_markers)\n",
    "method. Their appearance on the pattern can also be controlled to some extent.\n",
    "The above method itself calls\n",
    "[bands_as_markers()](../reference.rst#kikuchipy.simulations.GeometricalEBSDSimulation.bands_as_markers),\n",
    "[pc_as_markers()](../reference.rst#kikuchipy.simulations.GeometricalEBSDSimulation.pc_as_markers),\n",
    "[zone_axes_as_markers()](../reference.rst#kikuchipy.simulations.GeometricalEBSDSimulation.zone_axes_as_markers),\n",
    "and [zone_axes_labels_as_markers()](../reference.rst#kikuchipy.simulations.GeometricalEBSDSimulation.zone_axes_labels_as_markers).\n",
    "See their documentation for available modifications.\n",
    "\n",
    "Let's first remove the markers from the signal, and add only the bands and zone\n",
    "axes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del s.metadata.Markers\n",
    "s.add_marker(\n",
    "    marker=sim.as_markers(\n",
    "        pc=False,\n",
    "        zone_axes_labels=False,\n",
    "        bands_kwargs=dict(\n",
    "            family_colors=[\"w\", \"magenta\", \"cyan\", \"lime\"], linestyle=\"--\",\n",
    "        ),\n",
    "        zone_axes_kwargs=dict(\n",
    "            marker=\"s\", size=150, facecolor=\"none\", edgecolor=\"w\",\n",
    "        ),\n",
    "    ),\n",
    "    plot_marker=False,\n",
    "    permanent=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s.plot(navigator=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del s.metadata.Markers\n",
    "s.add_marker(marker=markers, plot_marker=False, permanent=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can write single EBSD patterns with the markers on top to file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = s._plot.signal_plot.figure\n",
    "bbox = matplotlib.transforms.Bbox.from_extents(\n",
    "    np.array(fig.axes[0].bbox.extents) / 72  # The denominator may vary\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nav_x, nav_y = s.axes_manager.indices\n",
    "temp_dir = tempfile.mkdtemp() + \"/\"\n",
    "fname = temp_dir + f\"geosim_y{nav_y}_x{nav_x}.png\"\n",
    "s._plot.signal_plot.figure.savefig(fname, bbox_inches=bbox, dpi=150)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbsphinx-thumbnail": {
     "tooltip": "Projection of the centers of Kikuchi bands and zone axes onto an EBSD detector"
    },
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [],
   "source": [
    "s.axes_manager.indices = (2, 1)\n",
    "s.plot(navigator=None, colorbar=False, axes_off=True, title=\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nav_x, nav_y = s.axes_manager.indices\n",
    "fname = temp_dir + f\"geosim_y{nav_y}_x{nav_x}.png\"\n",
    "s._plot.signal_plot.figure.savefig(fname, bbox_inches=bbox, dpi=150)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coordinates of these bands and zone axes are available as class attributes.\n",
    "For the bands, we can e.g. extract the plane trace coordinates (y0, x0, y1, x1)\n",
    "in either gnomonic or detector coordinates (taking into account the detector\n",
    "size in pixels) for all bands or per navigation position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.bands[0, 0].plane_trace_coordinates[:10]  # Gnomonic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.bands_detector_coordinates[0, 0, :10]  # Detector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The NaN values signify that that particular band is not visible on the detector\n",
    "in that position. The crystal plane normal of each band, pointing from the\n",
    "source point to the detector, is also available"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.bands[1, 1].hkl_detector[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And where the vector hits the detector, in either detector or gnomonic\n",
    "coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.bands[0, 0].x_detector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.bands[0, 0].x_gnomonic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same information is available for the zone axes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.zone_axes_label_detector_coordinates[0, 0][20:30]  # Detector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.zone_axes[0, 1].uvw_detector[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.zone_axes[0, 0].y_detector[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.zone_axes[0, 0].y_gnomonic[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this information, we can go around kikuchipy when plotting and\n",
    "only use matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nav_idx = (2, 1)[::-1]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "ax.imshow(s.inav[nav_idx], cmap=\"gray\")\n",
    "\n",
    "print(sim.bands_detector_coordinates.shape)\n",
    "for i in np.ndindex(sim.bands_detector_coordinates.shape[2]):\n",
    "    sim_slice = nav_idx + (i,)\n",
    "    coords = sim.bands_detector_coordinates[sim_slice][0]\n",
    "    y0, x0, y1, x1 = coords\n",
    "    ax.axline((y0, x0), (y1, x1), linestyle=\"--\", color=\"w\")\n",
    "\n",
    "print(sim.zone_axes_detector_coordinates.shape)\n",
    "for j in np.ndindex(sim.zone_axes_detector_coordinates.shape[2]):\n",
    "    sim_slice = nav_idx + (j,)\n",
    "    coords = sim.zone_axes_detector_coordinates[sim_slice][0]\n",
    "    x, y = coords\n",
    "    ax.scatter(x=x, y=y, zorder=5, s=50)\n",
    "\n",
    "_ = ax.axis((0, 59, 59, 0))\n",
    "ax.axis(\"off\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove files written to disk in this user guide\n",
    "import os\n",
    "for file in [\"geosim_y0_x0.png\", \"geosim_y1_x2.png\"]:\n",
    "    os.remove(temp_dir + file)\n",
    "os.rmdir(temp_dir)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}