{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cmp-title",
   "metadata": {},
   "source": [
    "# Comparing two REI point clouds\n",
    "\n",
    "`REIComparison` measures the spatial overlap of two rare-event point clouds\n",
    "(two metrics/thresholds, or a prediction against a reference), even when they sit\n",
    "on grids of different spacing (they must share an origin). It reports IoU, Dice,\n",
    "and containment, matches clusters one-to-one, and writes a classified VTK. Pure\n",
    "Python (numpy + scipy); runs from `pip install graintrace`.\n",
    "\n",
    "See :doc:`/algorithms/rei-comparison` and :class:`~graintrace.REIComparison`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cmp-imports",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-09T23:29:43.128663Z",
     "iopub.status.busy": "2026-10-09T23:29:43.127864Z",
     "iopub.status.idle": "2026-10-09T23:29:43.821106Z",
     "shell.execute_reply": "2026-10-09T23:29:43.819716Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from graintrace.rei_comparison import REIComparison\n",
    "\n",
    "import os\n",
    "\n",
    "os.makedirs(\"out/rei_cmp\", exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cmp-gen-md",
   "metadata": {},
   "source": [
    "## Two clouds of the same blobs at different spacing\n",
    "\n",
    "Both clouds sample the **same three spheres in physical units**; REI 1 uses a\n",
    "1 um lattice and REI 2 a 2 um lattice. Each sphere is one cluster id.\n",
    "\n",
    "This is the case the comparator is built for. The two grids have different\n",
    "spacing but share an origin, which is the only registration assumption\n",
    "`REIComparison` makes -- it applies no rotation and no translation, so clouds\n",
    "that are merely *near* each other will score as disagreeing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cmp-gen",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-09T23:29:43.825841Z",
     "iopub.status.busy": "2026-10-09T23:29:43.825625Z",
     "iopub.status.idle": "2026-10-09T23:29:43.844273Z",
     "shell.execute_reply": "2026-10-09T23:29:43.842908Z"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "os.makedirs(\"out/rei_cmp\", exist_ok=True)\n",
    "\n",
    "\n",
    "def blobs(extent, spacing, centers, radius):\n",
    "    \"\"\"Sample three spheres (physical units) onto a lattice of the given spacing.\"\"\"\n",
    "    g = np.arange(0.0, extent, spacing)\n",
    "    xx, yy, zz = np.meshgrid(g, g, g, indexing=\"ij\")\n",
    "    rows = []\n",
    "    for cid, (cx, cy, cz) in enumerate(centers, start=1):\n",
    "        m = (xx - cx) ** 2 + (yy - cy) ** 2 + (zz - cz) ** 2 <= radius * radius\n",
    "        rows.append(\n",
    "            pd.DataFrame(\n",
    "                {\"x\": xx[m], \"y\": yy[m], \"z\": zz[m], \"rare_cluster_id\": cid}\n",
    "            )\n",
    "        )\n",
    "    return pd.concat(rows, ignore_index=True)\n",
    "\n",
    "\n",
    "extent = 40.0\n",
    "centers = [(10.0, 10.0, 10.0), (25.0, 25.0, 20.0), (15.0, 30.0, 30.0)]\n",
    "df1 = blobs(extent, 1.0, centers, radius=6.0)   # fine grid, 1 um\n",
    "df2 = blobs(extent, 2.0, centers, radius=6.0)   # coarse grid, 2 um, same spheres\n",
    "df1.to_csv(\"out/rei_cmp/rei_A.csv\", index=False)\n",
    "df2.to_csv(\"out/rei_cmp/rei_B.csv\", index=False)\n",
    "print(\"REI 1 points:\", len(df1), \"| REI 2 points:\", len(df2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cmp-run-md",
   "metadata": {},
   "source": [
    "## Run the comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cmp-run",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-09T23:29:43.847770Z",
     "iopub.status.busy": "2026-10-09T23:29:43.847619Z",
     "iopub.status.idle": "2026-10-09T23:29:43.873648Z",
     "shell.execute_reply": "2026-10-09T23:29:43.872212Z"
    }
   },
   "outputs": [],
   "source": [
    "comp = REIComparison(\n",
    "    rei_csv_1=\"out/rei_cmp/rei_A.csv\", rei_csv_2=\"out/rei_cmp/rei_B.csv\",\n",
    "    output_dir=\"out/rei_cmp\",\n",
    "    spacing_1=1.0, spacing_2=2.0,\n",
    "    coord_cols=(\"x\", \"y\", \"z\"),\n",
    "    cluster_col=\"rare_cluster_id\",\n",
    ")\n",
    "result = comp.run_comparison()\n",
    "m = result[\"overlap_metrics\"] if \"overlap_metrics\" in result else result\n",
    "for k in (\"iou\", \"dice\", \"containment_1\", \"containment_2\"):\n",
    "    if k in m:\n",
    "        print(f\"{k:15} {m[k]:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ff9f4fe",
   "metadata": {},
   "source": [
    "### What the numbers mean\n",
    "\n",
    "All four metrics are computed on a common lattice: both clouds are resampled\n",
    "onto the finer spacing $s_{\\mathrm{ref}} = \\min(s_1, s_2)$ and the overlap is\n",
    "an exact integer-index intersection, so there is no tolerance and no\n",
    "nearest-neighbour matching involved.\n",
    "\n",
    "$$\\mathrm{IoU} = \\frac{|A \\cap B|}{|A \\cup B|},\n",
    "\\qquad\n",
    "\\mathrm{Dice} = \\frac{2|A \\cap B|}{|A| + |B|},\n",
    "\\qquad\n",
    "\\mathrm{containment}_1 = \\frac{|A \\cap B|}{|A|}.$$\n",
    "\n",
    "IoU and Dice are symmetric and answer \"do these agree?\". The two containments\n",
    "are not, and they are the pair that tells you *how* they disagree: a high\n",
    "$\\mathrm{containment}_1$ with a low $\\mathrm{containment}_2$ means region 1\n",
    "sits inside region 2, which is the signature of one criterion being strictly\n",
    "more permissive rather than simply wrong."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83192d0a",
   "metadata": {},
   "source": [
    "## Cluster-to-cluster matching\n",
    "\n",
    "Global overlap says whether the two regions coincide; it does not say whether\n",
    "they decomposed the region the same way. `cluster_match.csv` pairs clusters\n",
    "one-to-one by maximising total overlap volume (Hungarian assignment), so the\n",
    "pairing is label-agnostic \u2014 cluster 3 in one run can legitimately be cluster 1\n",
    "in the other. Unmatched clusters are flagged `-1`, which is how a split or a\n",
    "merge shows up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed48ca9d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-09T23:29:43.877405Z",
     "iopub.status.busy": "2026-10-09T23:29:43.876845Z",
     "iopub.status.idle": "2026-10-09T23:29:43.891679Z",
     "shell.execute_reply": "2026-10-09T23:29:43.890254Z"
    }
   },
   "outputs": [],
   "source": [
    "match = pd.read_csv(\"out/rei_cmp/cluster_match.csv\")\n",
    "print(match.to_string(index=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52f252f6",
   "metadata": {},
   "source": [
    "## Seeing the overlap\n",
    "\n",
    "`overlap_cloud.vtk` carries a `membership` scalar on the common lattice:\n",
    "1 = only in REI 1, 2 = only in REI 2, 3 = in both. ParaView is the natural\n",
    "viewer, but one slice through the classified cloud makes the point here.\n",
    "\n",
    "The disagreement is a **boundary** effect, and it goes both ways. The coarse\n",
    "cloud's voxels are 2 um cubes, so near the sphere surface each one either\n",
    "sticks out past the fine cloud's edge (only 2) or fails to reach it (only 1).\n",
    "Interior voxels agree exactly. That is the characteristic signature of a pure\n",
    "resolution difference, as opposed to two criteria flagging genuinely different\n",
    "places -- which would show up as whole clusters unmatched in the table above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f407e427",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-09T23:29:43.894978Z",
     "iopub.status.busy": "2026-10-09T23:29:43.894534Z",
     "iopub.status.idle": "2026-10-09T23:29:44.619077Z",
     "shell.execute_reply": "2026-10-09T23:29:44.617792Z"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import pyvista as pv\n",
    "\n",
    "cloud = pv.read(\"out/rei_cmp/overlap_cloud.vtk\")\n",
    "pts = np.asarray(cloud.points)\n",
    "mem = np.asarray(cloud[\"membership\"]).astype(int)\n",
    "\n",
    "zc = 10.0\n",
    "sel = pts[:, 2] == zc\n",
    "labels = {1: (\"only REI 1\", \"tab:blue\"), 2: (\"only REI 2\", \"tab:red\"), 3: (\"both\", \"tab:green\")}\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5.2, 5))\n",
    "for key, (lab, col) in labels.items():\n",
    "    m = sel & (mem == key)\n",
    "    ax.scatter(pts[m, 0], pts[m, 1], s=16, c=col, label=f\"{lab} ({int(m.sum())})\")\n",
    "ax.set_aspect(\"equal\")\n",
    "ax.set_xlabel(\"x [um]\")\n",
    "ax.set_ylabel(\"y [um]\")\n",
    "ax.set_title(f\"overlap_cloud.vtk membership, z = {zc:g} um slice\")\n",
    "ax.legend(loc=\"upper left\", fontsize=8)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cmp-seealso",
   "metadata": {},
   "source": [
    "The two clouds cover the same blobs, so IoU/Dice are high. `output_dir` also\n",
    "gets `overlap_metrics.json`, a classified `overlap_cloud.vtk` (only-1 / only-2 /\n",
    "both), and `cluster_match.csv` (the one-to-one cluster pairing).\n",
    "\n",
    "**See also**\n",
    "\n",
    "- Algorithm: :doc:`/algorithms/rei-comparison`\n",
    "- Producing REI point clouds: :doc:`rare-event-identification`"
   ]
  }
 ],
 "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}