{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial about computing hulls for localization data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each set of localizations with 2D or 3D spatial coordinates various hull can be computed. A hull can be the minimal bounding box, the oriented minimal bounding box, the convex hull, or an alpha shape. \n",
    "\n",
    "You can trigger computation of specific hull objects using a specific hull class or from the corresponding LocData attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import locan as lc\n",
    "from locan.data.hulls import BoundingBox, ConvexHull, OrientedBoundingBox"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lc.show_versions(system=False, dependencies=False, verbose=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Synthetic data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(seed=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "locdata = lc.simulate_Thomas(parent_intensity=1e-3, region=((0, 100), (0, 100)), cluster_mu=10, cluster_std=2, seed=rng)\n",
    "\n",
    "locdata.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1)\n",
    "locdata.data.plot.scatter(x='position_x', y='position_y', ax=ax, color='Blue', label='locdata')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Minimal bounding box for spatial coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# H = BoundingBox(locdata.coordinates)\n",
    "\n",
    "H = locdata.bounding_box\n",
    "\n",
    "print('dimension: ', H.dimension)\n",
    "print('hull: ', H.hull)\n",
    "print('width: ', H.width)\n",
    "print('vertices: ', H.vertices)\n",
    "print('region_measure: ', H.region_measure)\n",
    "print('subregion_measure: ', H.subregion_measure)\n",
    "print('region: ', H.region)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1)\n",
    "ax.add_patch(locdata.bounding_box.region.as_artist(alpha=0.2))\n",
    "locdata.data.plot.scatter(x='position_x', y='position_y', ax=ax, color='Blue', label='locdata')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Oriented minimal bounding box for spatial coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# H = OrientedBoundingBox(locdata.coordinates)\n",
    "\n",
    "H = locdata.oriented_bounding_box\n",
    "\n",
    "print('dimension: ', H.dimension)\n",
    "print('hull: ', H.hull)\n",
    "print('vertices: ', H.vertices)\n",
    "print('width: ', H.width)\n",
    "print('region_measure: ', H.region_measure)\n",
    "print('subregion_measure: ', H.subregion_measure)\n",
    "print('region: ', H.region)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1)\n",
    "ax.add_patch(locdata.oriented_bounding_box.region.as_artist(alpha=0.2))\n",
    "locdata.data.plot.scatter(x='position_x', y='position_y', ax=ax, color='Blue', label='locdata')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convex hull"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convex hull for spatial coordinates (scipy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# H = ConvexHull(locdata.coordinates, method='scipy')\n",
    "\n",
    "H = locdata.convex_hull\n",
    "\n",
    "print('dimension: ', H.dimension)\n",
    "print('hull: ', H.hull)\n",
    "# print('vertex_indices: ', H.vertex_indices)\n",
    "# print('vertices: ', H.vertices)\n",
    "print('region_measure: ', H.region_measure)\n",
    "print('subregion_measure: ', H.subregion_measure)\n",
    "print('points on boundary: ', H.points_on_boundary)\n",
    "print('points on boundary relative to all points: ', H.points_on_boundary_rel)\n",
    "print('region: ', H.region)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1)\n",
    "ax.add_patch(locdata.convex_hull.region.as_artist(alpha=0.2))\n",
    "locdata.data.plot.scatter(x='position_x', y='position_y', ax=ax, color='Blue', label='locdata')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convex hull for spatial coordinates (shapely)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some hulls can be computed from different algorithms. If implemented, use the `methods` parameter to specify the algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = ConvexHull(locdata.coordinates, method='shapely')\n",
    "\n",
    "print('dimension: ', H.dimension)\n",
    "print('hull: ', H.hull)\n",
    "# print('vertices: ', H.vertices)\n",
    "print('region_measure: ', H.region_measure)\n",
    "print('subregion_measure: ', H.subregion_measure)\n",
    "print('points on boundary: ', H.points_on_boundary)\n",
    "print('points on boundary relative to all points: ', H.points_on_boundary_rel)\n",
    "print('region: ', H.region)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1)\n",
    "ax.add_patch(H.region.as_artist(alpha=0.2))\n",
    "locdata.data.plot.scatter(x='position_x', y='position_y', ax=ax, color='Blue', label='locdata')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alpha shape for spatial coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The alpha shape depends on a single parameter `alpha` (not to confuse with the alpha to specify opacity in figures). The alpha complex is an alpha-independent representation of all alpha shapes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can get all apha values for which the corresponding alpha shape changes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lc.AlphaComplex(locdata.coordinates).alphas()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can determine an optimal `alpha`, i.e. the smallest `alpha` for which all points are still part of the alpha shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt_alpha = lc.AlphaComplex(locdata.coordinates).optimal_alpha()\n",
    "opt_alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# H = lc.AlphaShape(opt_alpha, locdata.coordinates)\n",
    "\n",
    "locdata.update_alpha_shape(alpha=opt_alpha)\n",
    "H = locdata.alpha_shape\n",
    "\n",
    "print('dimension: ', H.dimension)\n",
    "# print('vertex_indices: ', H.vertex_indices)\n",
    "# print('vertices: ', H.vertices)\n",
    "print('region_measure: ', H.region_measure)\n",
    "print('subregion_measure: ', H.subregion_measure)\n",
    "print('points in alpha shape: ', H.n_points_alpha_shape)\n",
    "print('points in alpha shape relative to all points: ', H.n_points_alpha_shape_rel)\n",
    "print('points on boundary: ', H.n_points_on_boundary)\n",
    "print('points on boundary relative to all points: ', H.n_points_on_boundary_rel)\n",
    "print('region: ', H.region)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The alpha shape is made of different vertex types that can be differentiated as *exterior*, *interior*, *regular* or *singular*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ac_simplices_all = H.alpha_complex.get_alpha_complex_lines(H.alpha, type='all')\n",
    "ac_simplices_exterior = H.alpha_complex.get_alpha_complex_lines(H.alpha, type='exterior')\n",
    "ac_simplices_interior = H.alpha_complex.get_alpha_complex_lines(H.alpha, type='interior')\n",
    "ac_simplices_regular = H.alpha_complex.get_alpha_complex_lines(H.alpha, type='regular')\n",
    "ac_simplices_singular = H.alpha_complex.get_alpha_complex_lines(H.alpha, type='singular')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1)\n",
    "\n",
    "for simp in ac_simplices_all:\n",
    "    ax.plot(locdata.coordinates[simp, 0], locdata.coordinates[simp, 1], '-b')\n",
    "for simp in ac_simplices_interior:\n",
    "    ax.plot(locdata.coordinates[simp, 0], locdata.coordinates[simp, 1], '--g')\n",
    "for simp in ac_simplices_regular:\n",
    "    ax.plot(locdata.coordinates[simp, 0], locdata.coordinates[simp, 1], '--r')\n",
    "for simp in ac_simplices_singular:\n",
    "    ax.plot(locdata.coordinates[simp, 0], locdata.coordinates[simp, 1], '--y')\n",
    "\n",
    "locdata.data.plot.scatter(x='position_x', y='position_y', ax=ax, color='Blue', label='locdata')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Often the *regular* representation is good enough."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1)\n",
    "for simp in ac_simplices_regular:\n",
    "    ax.plot(locdata.coordinates[simp, 0], locdata.coordinates[simp, 1], '-r')\n",
    "locdata.data.plot.scatter(x='position_x', y='position_y', ax=ax, color='Blue', label='locdata')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can get the connected components as list of `Region`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H.connected_components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "connected_component_0 = H.connected_components[0]\n",
    "\n",
    "print('dimension: ', connected_component_0.dimension)\n",
    "print('region_measure: ', connected_component_0.region_measure)\n",
    "print('subregion_measure: ', connected_component_0.subregion_measure)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The alpha shape for a smaller alpha can have multiple connected components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = lc.AlphaShape(5, locdata.coordinates)\n",
    "\n",
    "print('dimension: ', H.dimension)\n",
    "# print('vertex_indices: ', H.vertex_indices)\n",
    "# print('vertices: ', H.vertices)\n",
    "print('region_measure: ', H.region_measure)\n",
    "print('subregion_measure: ', H.subregion_measure)\n",
    "print('points in alpha shape: ', H.n_points_alpha_shape)\n",
    "print('points in alpha shape relative to all points: ', H.n_points_alpha_shape_rel)\n",
    "print('points on boundary: ', H.n_points_on_boundary)\n",
    "print('points on boundary relative to all points: ', H.n_points_on_boundary_rel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ac_simplices_all = H.alpha_complex.get_alpha_complex_lines(H.alpha, type='all')\n",
    "ac_simplices_exterior = H.alpha_complex.get_alpha_complex_lines(H.alpha, type='exterior')\n",
    "ac_simplices_interior = H.alpha_complex.get_alpha_complex_lines(H.alpha, type='interior')\n",
    "ac_simplices_regular = H.alpha_complex.get_alpha_complex_lines(H.alpha, type='regular')\n",
    "ac_simplices_singular = H.alpha_complex.get_alpha_complex_lines(H.alpha, type='singular')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1)\n",
    "\n",
    "for simp in ac_simplices_all:\n",
    "    ax.plot(locdata.coordinates[simp, 0], locdata.coordinates[simp, 1], '-b')\n",
    "for simp in ac_simplices_interior:\n",
    "    ax.plot(locdata.coordinates[simp, 0], locdata.coordinates[simp, 1], '--g')\n",
    "for simp in ac_simplices_regular:\n",
    "    ax.plot(locdata.coordinates[simp, 0], locdata.coordinates[simp, 1], '--r')\n",
    "for simp in ac_simplices_singular:\n",
    "    ax.plot(locdata.coordinates[simp, 0], locdata.coordinates[simp, 1], '--y')\n",
    "\n",
    "locdata.data.plot.scatter(x='position_x', y='position_y', ax=ax, color='Blue', label='locdata')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *regular* representation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1)\n",
    "for simp in ac_simplices_regular:\n",
    "    ax.plot(locdata.coordinates[simp, 0], locdata.coordinates[simp, 1], '-r')\n",
    "locdata.data.plot.scatter(x='position_x', y='position_y', ax=ax, color='Blue', label='locdata')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The connected components:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H.connected_components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "connected_component_0 = H.connected_components[0]\n",
    "\n",
    "print('dimension: ', connected_component_0.dimension)\n",
    "print('region_measure: ', connected_component_0.region_measure)\n",
    "print('subregion_measure: ', connected_component_0.subregion_measure)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.8"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}