{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden",
    "tags": []
   },
   "source": [
    "This notebook is part of the *orix* documentation https://orix.readthedocs.io. Links to the documentation won’t work from the notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Clustering across fundamental region boundaries\n",
    "\n",
    "In this tutorial we will perform density based clustering of crystal orientations with and without the application of crystal symmetry symmetry using simulated data, as presented in <cite data-cite=\"johnstone2020density\">Johnstone et al. (2020)</cite>.\n",
    "\n",
    "Import orix classes and various dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Exchange inline for notebook (or qt5 from pyqt) for interactive plotting\n",
    "%matplotlib inline\n",
    "\n",
    "# Import core external\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.cluster import DBSCAN\n",
    "\n",
    "# Colorisation & Animation\n",
    "from skimage.color import label2rgb\n",
    "from matplotlib.colors import to_rgb\n",
    "import matplotlib.animation as animation\n",
    "\n",
    "# Import orix classes\n",
    "from orix.quaternion import Orientation, OrientationRegion, Rotation\n",
    "from orix.quaternion.symmetry import C1, Oh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Generate artificial data\n",
    "\n",
    "Generate three random von Mises distributions of orientations as model clusters and set the *Oh* ($m\\bar{3}m$) point group symmetry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "n_orientations = 50\n",
    "alpha = 65  # Lower value gives \"looser\" distribution\n",
    "\n",
    "# Cluster 1\n",
    "cluster1 = Rotation.random_vonmises(n_orientations, alpha=alpha)\n",
    "\n",
    "# Cluster 2\n",
    "centre2 = Rotation.from_axes_angles((1, 0, 0), np.pi / 4)\n",
    "cluster2 = Rotation.random_vonmises(\n",
    "    n_orientations, alpha=alpha, reference=centre2\n",
    ")\n",
    "\n",
    "# Cluster 3\n",
    "centre3 = Rotation.from_axes_angles((1, 1, 0), np.pi / 3)\n",
    "cluster3 = Rotation.random_vonmises(\n",
    "    n_orientations, alpha=alpha, reference=centre3\n",
    ")\n",
    "\n",
    "# Stack and map into the Oh fundamental zone\n",
    "ori = Orientation.stack([cluster1, cluster2, cluster3]).flatten()\n",
    "ori.symmetry = Oh\n",
    "ori = ori.map_into_symmetry_reduced_zone()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Orientation clustering\n",
    "\n",
    "### Perform clustering without application of crystal symmetry\n",
    "\n",
    "Compute misorientations, i.e. distance between orientations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Remove symmetry by setting it to point group 1 (identity operation)\n",
    "ori_without_symmetry = Orientation(ori.data, symmetry=C1)\n",
    "\n",
    "# Misorientations\n",
    "mori1 = (~ori_without_symmetry).outer(ori_without_symmetry)\n",
    "\n",
    "# Misorientation angles\n",
    "D1 = mori1.angle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dbscan_naive = DBSCAN(eps=0.3, min_samples=10, metric=\"precomputed\").fit(D1)\n",
    "print(\"Labels:\", np.unique(dbscan_naive.labels_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perform clustering with application of crystal symmetry\n",
    "\n",
    "Compute misorientations, i.e. distance between orientations, with symmetry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mori2 = (~ori).outer(ori)\n",
    "\n",
    "mori2.symmetry = Oh\n",
    "mori2 = mori2.map_into_symmetry_reduced_zone()\n",
    "\n",
    "D2 = mori2.angle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dbscan = DBSCAN(\n",
    "    eps=np.deg2rad(17), min_samples=20, metric=\"precomputed\"\n",
    ").fit(D2.astype(np.float32))\n",
    "print(\"Labels:\", np.unique(dbscan.labels_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This should have shown that without symmetry there are 6 clusters, whereas with symmetry there are 3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualisation\n",
    "\n",
    "Assign colours to each cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_names = [to_rgb(f\"C{i}\") for i in range(6)]  # ['C0', 'C1', ...]\n",
    "\n",
    "colors_naive = label2rgb(\n",
    "    dbscan_naive.labels_, colors=color_names, bg_label=-1\n",
    ")\n",
    "colors = label2rgb(dbscan.labels_, colors=color_names, bg_label=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot orientation clusters with Matplotlib and [(Mis)orientation.scatter()](../reference/generated/orix.quaternion.Misorientation.scatter.rst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbsphinx-thumbnail": {
     "tooltip": "Density based clustering of crystal orientations with and without crystal symmetry"
    },
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [],
   "source": [
    "# Set symmetry to \"trick\" the scatter plot to use the Oh fundamental zone\n",
    "ori_without_symmetry.symmetry = ori.symmetry\n",
    "\n",
    "# Create figure with a height/width ratio of 1/2\n",
    "fig = plt.figure(figsize=(12, 6))\n",
    "\n",
    "# Add the fundamental zones with clusters to the existing figure\n",
    "ori_without_symmetry.scatter(figure=fig, position=(1, 2, 1), c=colors_naive)\n",
    "ori.scatter(figure=fig, position=122, c=colors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate an animation of the plot (assuming an interactive Matplotlib backend is used)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def animate(angle):\n",
    "    fig.axes[0].view_init(15, angle)\n",
    "    fig.axes[1].view_init(15, angle)\n",
    "    plt.draw()\n",
    "\n",
    "\n",
    "ani = animation.FuncAnimation(\n",
    "    fig, animate, np.linspace(75, 360 + 74, 720), interval=25\n",
    ")"
   ]
  }
 ],
 "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}