{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial about radial distance distribution analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pairwise distance distribution p(r) - as derived from a histogram of\n",
    "pairwise distances - represents the probability distribution function\n",
    "to find for a localization at r = 0 another localization at distance r + delta_r.\n",
    "\n",
    "The radial distribution function (also called pair correlation function) g(r)\n",
    "represents the pairwise distance distribution function p(r) normalized to the\n",
    "region measure of the circular ring or shell with inner radius r and outer\n",
    "radius r + delta_r and relative to the expected number of localizations per unit area for\n",
    "complete spatial randomness."
   ]
  },
  {
   "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 matplotlib.pyplot as plt\n",
    "\n",
    "import locan as lc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lc.show_versions(system=False, dependencies=False, verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(seed=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Synthetic data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We simulate localization data at two different intensities (localization density) that is (i) homogeneously Poisson distributed (also described as complete spatial randomness, csr) and that (ii) follows a Neyman-Scott distribution (blobs)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "locdata_csr_0 = lc.simulate_Poisson(intensity=1e-3, region=((0,1000), (0,1000)), seed=rng)\n",
    "locdata_csr_1 = lc.simulate_Poisson(intensity=1e-2, region=((0,1000), (0,1000)), seed=rng)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "locdata_blob_0 = lc.simulate_Thomas(parent_intensity=1e-4, region=((0, 1000), (0, 1000)), cluster_mu=10, cluster_std=5, seed=rng)\n",
    "locdata_blob_1 = lc.simulate_Thomas(parent_intensity=1e-3, region=((0, 1000), (0, 1000)), cluster_mu=10, cluster_std=5, seed=rng)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Number of localizations:\")\n",
    "print(\"csr_0:\", len(locdata_csr_0))\n",
    "print(\"csr_1:\", len(locdata_csr_1))\n",
    "print(\"blob_0:\", len(locdata_blob_0))\n",
    "print(\"blob_1:\", len(locdata_blob_1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scatter plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=2, ncols=2)\n",
    "locdata_csr_0.data.plot.scatter(x='position_x', y='position_y', ax=axes[0, 0], color='Blue', s=1, alpha=0.1, label='locdata_csr')\n",
    "locdata_csr_1.data.plot.scatter(x='position_x', y='position_y', ax=axes[0, 1], color='Blue', s=1, alpha=0.1, label='locdata_csr')\n",
    "locdata_blob_0.data.plot.scatter(x='position_x', y='position_y', ax=axes[1, 0], color='Blue', s=1, alpha=0.1, label='locdata_blobs')\n",
    "locdata_blob_1.data.plot.scatter(x='position_x', y='position_y', ax=axes[1, 1], color='Blue', s=1, alpha=0.1, label='locdata_blobs')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Radial distance distribution "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We determine all pairwise distances and the normalized radial distance distribution (resembling the pair correlation function)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = np.linspace(0, 300, 100)\n",
    "rdf_csr_0 = lc.RadialDistribution(bins=bins).compute(locdata_csr_0)\n",
    "rdf_csr_1 = lc.RadialDistribution(bins=bins).compute(locdata_csr_1)\n",
    "rdf_blob_0 = lc.RadialDistribution(bins=bins).compute(locdata_blob_0)\n",
    "rdf_blob_1 = lc.RadialDistribution(bins=bins).compute(locdata_blob_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rdf_csr_0.hist(alpha=0.5, label=\"csr_0\")\n",
    "rdf_blob_0.hist(alpha=0.5, label=\"blob_0\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rdf_csr_1.hist(alpha=0.5, label=\"csr_1\")\n",
    "rdf_blob_1.hist(alpha=0.5, label=\"blob_1\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Synthetic data - 3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "locdata_csr_0 = lc.simulate_Poisson(intensity=1e-6, region=((0,1000), (0,1000), (0,1000)), seed=rng)\n",
    "locdata_csr_1 = lc.simulate_Poisson(intensity=1e-5, region=((0,1000), (0,1000), (0,1000)), seed=rng)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "locdata_blob_0 = lc.simulate_Thomas(parent_intensity=1e-7, region=((0, 1000), (0, 1000), (0,1000)), cluster_mu=10, cluster_std=5, seed=rng)\n",
    "locdata_blob_1 = lc.simulate_Thomas(parent_intensity=1e-6, region=((0, 1000), (0, 1000), (0,1000)), cluster_mu=10, cluster_std=5, seed=rng)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Number of localizations:\")\n",
    "print(\"csr_0:\", len(locdata_csr_0))\n",
    "print(\"csr_1:\", len(locdata_csr_1))\n",
    "print(\"blob_0:\", len(locdata_blob_0))\n",
    "print(\"blob_1:\", len(locdata_blob_1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We determine all pairwise distances and normalized radial distance distribution (resembling the pair correlation function)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = np.linspace(0, 300, 100)\n",
    "rdf_csr_0 = lc.RadialDistribution(bins=bins).compute(locdata_csr_0)\n",
    "rdf_csr_1 = lc.RadialDistribution(bins=bins).compute(locdata_csr_1)\n",
    "rdf_blob_0 = lc.RadialDistribution(bins=bins).compute(locdata_blob_0)\n",
    "rdf_blob_1 = lc.RadialDistribution(bins=bins).compute(locdata_blob_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rdf_csr_0.hist(alpha=0.5, label=\"csr_0\")\n",
    "rdf_csr_1.hist(alpha=0.5, label=\"csr_1\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rdf_csr_0.hist(alpha=0.5, label=\"csr_0\")\n",
    "rdf_blob_0.hist(alpha=0.5, label=\"blob_0\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rdf_csr_1.hist(alpha=0.5, label=\"csr_1\")\n",
    "rdf_blob_1.hist(alpha=0.5, label=\"blob_1\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Radial distance distribution with precomputed distances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use precomputed pairwise distances for quick variation of bin (radii) settings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "locdata_blob_3D = lc.simulate_Thomas(parent_intensity=1e-7, region=((0, 1000), (0, 1000), (0, 1000)), cluster_mu=10, cluster_std=5, seed=rng)\n",
    "locdata_blob_3D.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd_blob_3D = lc.PairDistances().compute(locdata_blob_3D)\n",
    "pd_blob_3D.results.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = np.linspace(0, 100, 25)\n",
    "rdf_blob_3D_0 = lc.RadialDistribution(bins=bins, pair_distances=pd_blob_3D).compute(locdata_blob_3D)\n",
    "bins = np.linspace(0, 100, 100)\n",
    "rdf_blob_3D_1 = lc.RadialDistribution(bins=bins, pair_distances=pd_blob_3D).compute(locdata_blob_3D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rdf_blob_3D_0.hist(alpha=0.5, label=\"blob_3D_0\")\n",
    "rdf_blob_3D_1.hist(alpha=0.5, label=\"blob_3D_1\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Radial distribution function for various regions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compare radial distribution functions for localizations homogeneously distributed in bounded regions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = np.linspace(0, 1000, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region_rectangle = lc.Rectangle((0, 0), 1000, 1000, 30)\n",
    "locdata_rectangle = lc.simulate_Poisson(intensity=1e-3, region=region_rectangle, seed=rng)\n",
    "rdf_rectangle = lc.RadialDistribution(bins=bins).compute(locdata=locdata_rectangle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1)\n",
    "locdata_rectangle.data.plot.scatter(x='position_x', y='position_y', ax=ax, color='Blue', label='locdata')\n",
    "locdata_rectangle.region.plot(ax=ax, fill=False, color='Red')\n",
    "ax.set_aspect(\"equal\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region_circle = lc.Ellipse((500, 500), 1000, 1000, 0)\n",
    "locdata_circle = lc.simulate_Poisson(intensity=1e-3, region=region_circle, seed=rng)\n",
    "rdf_circle = lc.RadialDistribution(bins=bins).compute(locdata=locdata_circle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1)\n",
    "locdata_circle.data.plot.scatter(x='position_x', y='position_y', ax=ax, color='Blue', label='locdata')\n",
    "locdata_circle.region.plot(ax=ax, fill=False, color='Red')\n",
    "ax.set_aspect(\"equal\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region_ellipse = lc.Ellipse((500, 500), 100, 1000, 45)\n",
    "locdata_ellipse = lc.simulate_Poisson(intensity=1e-2, region=region_ellipse, seed=rng)\n",
    "rdf_ellipse = lc.RadialDistribution(bins=bins).compute(locdata=locdata_ellipse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1)\n",
    "locdata_ellipse.data.plot.scatter(x='position_x', y='position_y', ax=ax, color='Blue', label='locdata')\n",
    "locdata_ellipse.region.plot(ax=ax, fill=False, color='Red');\n",
    "ax.set_aspect(\"equal\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region_cuboid = lc.AxisOrientedCuboid((0, 0, 0), 1000, 1000, 1000)\n",
    "locdata_cuboid = lc.simulate_Poisson(intensity=1e-6, region=region_cuboid, seed=rng)\n",
    "rdf_cuboid = lc.RadialDistribution(bins=bins).compute(locdata=locdata_cuboid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1)\n",
    "locdata_cuboid.data.plot.scatter(x='position_x', y='position_y', ax=ax, color='Blue', label='locdata')\n",
    "#locdata_cuboid.region.plot(ax=ax, fill=False, color='Red');\n",
    "ax.set_aspect(\"equal\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rdf_rectangle.hist(label=\"rectangle\", histtype=\"step\")\n",
    "rdf_circle.hist(label=\"circle\", histtype=\"step\")\n",
    "rdf_ellipse.hist(label=\"ellipse\", histtype=\"step\")\n",
    "rdf_cuboid.hist(label=\"cuboid\", histtype=\"step\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Radial distribution function normalized for specific regions of interest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To correct for boundary effects from the enclosing region, compute a radial distance distribution relative to that of a homogeneous point process in the same region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region_ellipse = lc.Ellipse((500, 500), 300, 1000, 45)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "locdata_csr = lc.simulate_Poisson(intensity=1e-2, region=region_ellipse, seed=rng)\n",
    "locdata_blobs = lc.simulate_Thomas(parent_intensity=1e-3, region=region_ellipse, cluster_mu=10, cluster_std=5, seed=rng)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=2)\n",
    "locdata_csr.data.plot.scatter(x='position_x', y='position_y', ax=axes[0], color='Blue', s=1, alpha=0.1, label='locdata_csr')\n",
    "locdata_csr.region.plot(ax=axes[0], fill=False, color='Red')\n",
    "locdata_blobs.data.plot.scatter(x='position_x', y='position_y', ax=axes[1], color='Blue', s=1, alpha=0.1, label='locdata_blobs')\n",
    "locdata_blobs.region.plot(ax=axes[1], fill=False, color='Red')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "pd_csr = lc.PairDistances().compute(locdata=locdata_csr)\n",
    "pd_blob = lc.PairDistances().compute(locdata=locdata_blob)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = np.linspace(0, 1000, 100)\n",
    "rdf_csr = lc.RadialDistribution(bins=bins).compute(locdata=locdata_csr)\n",
    "rdf_blobs = lc.RadialDistribution(bins=bins).compute(locdata=locdata_blobs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rdf_csr.hist(alpha=0.5, label=\"csr\")\n",
    "rdf_blobs.hist(alpha=0.5, label=\"blobs\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a new RadialDistribution instance to carry the relative data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_corrected = lc.RadialDistributionResults()\n",
    "results_corrected.radii = rdf_blobs.results.radii\n",
    "results_corrected.data = rdf_blobs.results.data / rdf_csr.results.data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rdf = lc.RadialDistribution()\n",
    "rdf.results = results_corrected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rdf.hist(label=\"blobs vs csr\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The peak at small radii represents the cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = np.linspace(0, 100, 100)\n",
    "rdf_csr = lc.RadialDistribution(bins=bins).compute(locdata=locdata_csr)\n",
    "rdf_blobs = lc.RadialDistribution(bins=bins).compute(locdata=locdata_blobs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rdf_csr.hist(alpha=0.5, label=\"csr\")\n",
    "rdf_blobs.hist(alpha=0.5, label=\"blobs\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.12.6"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}