{
 "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": {},
   "source": [
    "# Clustering misorientations\n",
    "\n",
    "In this tutorial we will cluster Ti crystal misorientations using data obtained from a highly deformed specimen, using EBSD, as presented in <cite data-cite=\"johnstone2020density\">Johnstone et al. (2020)</cite>. The data can be downloaded to your local cache via the [orix.data](../reference/generated/orix.data.rst) module.\n",
    "\n",
    "Import orix classes and various dependencies"
   ]
  },
  {
   "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",
    "from matplotlib.colors import to_rgb\n",
    "from matplotlib.lines import Line2D\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from skimage.color import label2rgb\n",
    "from sklearn.cluster import DBSCAN\n",
    "\n",
    "from orix import data, plot\n",
    "from orix.quaternion import Orientation, Misorientation, Rotation\n",
    "from orix.quaternion.symmetry import D6\n",
    "from orix.vector import Vector3d\n",
    "\n",
    "\n",
    "plt.rcParams.update({\"font.size\": 20, \"figure.figsize\": (10, 10)})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import data\n",
    "\n",
    "Load Ti orientations with the point group symmetry *D6* (*622*). We have to explicitly allow downloading from an external source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ori = data.ti_orientations(allow_download=True)\n",
    "ori"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The orientations define transformations from the sample (lab) to the crystal reference frame, i.e. the Bunge convention. The above referenced paper assumes the opposite convention, which is the one used in MTEX. So, we have to invert the orientations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ori = ~ori"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reshape the orientation mapping data to the correct spatial dimension for the scan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ori = ori.reshape(381, 507)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select a subset of the orientations with a suitable size for this demonstration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ori = ori[-100:, :200]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot orientation maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ckey = plot.IPFColorKeyTSL(D6)\n",
    "\n",
    "directions = [(1, 0, 0), (0, 1, 0)]\n",
    "titles = [\"X\", \"Y\"]\n",
    "\n",
    "fig, axes = plt.subplots(ncols=2, figsize=(15, 10))\n",
    "for i, ax in enumerate(axes):\n",
    "    ckey.direction = Vector3d(directions[i])\n",
    "    # Invert because orix assumes lab2crystal when coloring orientations\n",
    "    ax.imshow(ckey.orientation2color(~ori))\n",
    "    ax.set_title(f\"IPF-{titles[i]}\")\n",
    "    ax.axis(\"off\")\n",
    "\n",
    "# Add color key\n",
    "ax_ipfkey = fig.add_axes(\n",
    "    [0.932, 0.37, 0.1, 0.1],  # (Left, bottom, width, height)\n",
    "    projection=\"ipf\",\n",
    "    symmetry=ori.symmetry.laue,\n",
    ")\n",
    "ax_ipfkey.plot_ipf_color_key()\n",
    "ax_ipfkey.set_title(\"\")\n",
    "fig.subplots_adjust(wspace=0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Map the orientations into the fundamental zone (find symmetrically equivalent orientations with the smallest angle of rotation) of *D6*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ori = ori.map_into_symmetry_reduced_zone()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute misorientations (in the horizontal direction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mori_all = Misorientation(~ori[:, :-1] * ori[:, 1:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Keep only misorientations with a disorientation angle higher than 7$^{\\circ}$, assumed to represent grain boundaries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boundary_mask = mori_all.angle > np.deg2rad(7)\n",
    "mori = mori_all[boundary_mask]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Map the misorientations into the fundamental zone of (*D6*, *D6*)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mori.symmetry = (D6, D6)\n",
    "mori = mori.map_into_symmetry_reduced_zone()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute distance matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Increase the chunk size for a faster but more memory intensive computation\n",
    "D = mori.get_distance_matrix()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clustering\n",
    "\n",
    "Apply mask to remove small misorientations associated with grain orientation spread"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "small_mask = mori.angle < np.deg2rad(7)\n",
    "D = D[~small_mask][:, ~small_mask]\n",
    "mori = mori[~small_mask]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For parameter explanations of the DBSCAN algorithm (Density-Based Spatial Clustering for Applications with Noise), see the [scikit-learn documentation](https://scikit-learn.org/stable/modules/generated/sklearn.cluster.DBSCAN.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute clusters\n",
    "dbscan = DBSCAN(eps=0.05, min_samples=10, metric=\"precomputed\").fit(D)\n",
    "\n",
    "unique_labels, all_cluster_sizes = np.unique(\n",
    "    dbscan.labels_, return_counts=True\n",
    ")\n",
    "print(\"Labels:\", unique_labels)\n",
    "\n",
    "n_clusters = unique_labels.size - 1\n",
    "print(\"Number of clusters:\", n_clusters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the mean misorientation associated with each cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_cluster_labels = unique_labels[\n",
    "    1:\n",
    "]  # Without the \"no-cluster\" label -1\n",
    "cluster_sizes = all_cluster_sizes[1:]\n",
    "\n",
    "rc = Rotation.from_axes_angles((0, 0, 1), 15, degrees=True)\n",
    "\n",
    "mori_mean = []\n",
    "for label in unique_cluster_labels:\n",
    "    # Rotate\n",
    "    mori_i = rc * mori[dbscan.labels_ == label]\n",
    "\n",
    "    # Map into the fundamental zone\n",
    "    mori_i.symmetry = (D6, D6)\n",
    "    mori_i = mori_i.map_into_symmetry_reduced_zone()\n",
    "\n",
    "    # Get the cluster mean\n",
    "    mori_i = mori_i.mean()\n",
    "\n",
    "    # Rotate back and add to list\n",
    "    cluster_mean_local = (~rc) * mori_i\n",
    "    mori_mean.append(cluster_mean_local)\n",
    "\n",
    "cluster_means = Misorientation.stack(mori_mean).flatten()\n",
    "\n",
    "# Map into the fundamental zone\n",
    "cluster_means.symmetry = (D6, D6)\n",
    "cluster_means = cluster_means.map_into_symmetry_reduced_zone()\n",
    "cluster_means"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspect misorientations in the axis-angle representation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_means.axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.rad2deg(cluster_means.angle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define reference misorientations associated with twinning orientation relationships"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# From Krakow et al.\n",
    "twin_theory = Rotation.from_axes_angles(\n",
    "    axes=[\n",
    "        (1, 0, 0),  # sigma7a\n",
    "        (1, 0, 0),  # sigma11a\n",
    "        (2, 1, 0),  # sigma11b\n",
    "        (1, 0, 0),  # sigma13a\n",
    "        (2, 1, 0),  # sigma13b\n",
    "    ],\n",
    "    angles=[64.40, 34.96, 85.03, 76.89, 57.22],\n",
    "    degrees=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate difference, defined as minimum rotation angle, between measured and theoretical values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mori2 = (~twin_theory).outer(cluster_means)\n",
    "sym_op = D6.outer(D6).unique()\n",
    "mori2_equiv = (\n",
    "    D6.outer(~twin_theory).outer(sym_op).outer(cluster_means).outer(D6)\n",
    ")\n",
    "D2 = mori2_equiv.angle.min(axis=(0, 2, 4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.rad2deg(D2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the first, second, third, and fourth clusters are within $4.5^{\\circ}$ of $\\Sigma7$a, $\\Sigma13$a, $\\Sigma11$a, and $\\Sigma13$b, respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualisation\n",
    "\n",
    "Associate colours with clusters for plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = [to_rgb(f\"C{i}\") for i in range(cluster_means.size)]\n",
    "labels_rgb = label2rgb(dbscan.labels_, colors=colors, bg_label=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspect misorientation axes of clusters in an inverse pole figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_sizes = all_cluster_sizes[1:]\n",
    "cluster_sizes_scaled = 1000 * cluster_sizes / cluster_sizes.max()\n",
    "\n",
    "fig, ax = plt.subplots(\n",
    "    figsize=(5, 5), subplot_kw=dict(projection=\"ipf\", symmetry=D6)\n",
    ")\n",
    "ax.scatter(\n",
    "    cluster_means.axis, c=colors, s=cluster_sizes_scaled, alpha=0.5, ec=\"k\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot a top view of the misorientation clusters within the fundamental zone for the (*D6*, *D6*) bicrystal symmetry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbsphinx-thumbnail": {
     "tooltip": "Density based clustering of crystal misorientations"
    },
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [],
   "source": [
    "wireframe_kwargs = dict(\n",
    "    color=\"black\", linewidth=0.5, alpha=0.1, rcount=361, ccount=361\n",
    ")\n",
    "fig = mori.scatter(\n",
    "    projection=\"axangle\",\n",
    "    wireframe_kwargs=wireframe_kwargs,\n",
    "    c=labels_rgb.reshape(-1, 3),\n",
    "    s=4,\n",
    "    alpha=0.5,\n",
    "    return_figure=True,\n",
    ")\n",
    "ax = fig.axes[0]\n",
    "ax.view_init(elev=90, azim=-60)\n",
    "\n",
    "handle_kwds = dict(marker=\"o\", color=\"none\", markersize=10)\n",
    "handles = []\n",
    "for i in range(n_clusters):\n",
    "    line = Line2D(\n",
    "        [0], [0], label=i + 1, markerfacecolor=colors[i], **handle_kwds\n",
    "    )\n",
    "    handles.append(line)\n",
    "ax.legend(\n",
    "    handles=handles,\n",
    "    loc=\"upper left\",\n",
    "    numpoints=1,\n",
    "    labelspacing=0.15,\n",
    "    columnspacing=0.15,\n",
    "    handletextpad=0.05,\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot side view of misorientation clusters in the fundamental zone for the (*D6*, *D6*) bicrystal symmetry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig2 = mori.scatter(\n",
    "    return_figure=True,\n",
    "    wireframe_kwargs=wireframe_kwargs,\n",
    "    c=labels_rgb.reshape(-1, 3),\n",
    "    s=4,\n",
    ")\n",
    "ax2 = fig2.axes[0]\n",
    "ax2.view_init(elev=0, azim=-60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate map of boundaries colored according to cluster membership"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mapping = np.ones(mori_all.shape + (3,))\n",
    "new_mask = (\n",
    "    np.where(boundary_mask)[0][~small_mask],\n",
    "    np.where(boundary_mask)[1][~small_mask],\n",
    ")\n",
    "mapping[new_mask] = labels_rgb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot map of boundaries colored according to cluster membership"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig3, ax3 = plt.subplots(figsize=(15, 10))\n",
    "ax3.imshow(mapping)\n",
    "ax3.set_xticks([])\n",
    "ax3.set_yticks([]);"
   ]
  }
 ],
 "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
}