exercise1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a39217e9-0d31-488c-9be3-8edc3315717a",
   "metadata": {},
   "source": [
    "# Exercise 1/3: Tracking by detection and simple frame-by-frame matching\n",
    "\n",
    "You could also run this notebook on your laptop, a GPU is not needed.\n",
    "\n",
    "Here we will walk through all basic components of a tracking-by-detection algorithm.\n",
    "\n",
    "You will learn\n",
    "- to **store and visualize** tracking results with `napari` (Exercise 1.1).\n",
    "- to use a robust pretrained deep-learning-based **object detection** algorithm called *StarDist* (Exercise 1.2).\n",
    "- to implement a basic **nearest-neighbor linking algorithm** (Exercises 1.3 - 1.6).\n",
    "- to compute optimal frame-by-frame linking by setting up a **bipartite matching problem** and using a python-based solver (Exercise 1.7).\n",
    "- to compute suitable object **features** for the object linking process with `scikit-image` (Exercise 1.8).\n",
    "\n",
    "Places where you are expected to write code are marked with ```YOUR CODE HERE```."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "014c728f-b294-4d20-83a9-48acbd227b62",
   "metadata": {},
   "source": [
    "![SegmentLocal](figures/trackmate-stardist-tracking.gif \"segment\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fee04091-df5a-43f2-bed4-a8643b44127b",
   "metadata": {
    "incorrectly_encoded_metadata": "tags=[] jp-MarkdownHeadingCollapsed=true",
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22bdc7a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Force keras to run on CPU\n",
    "import os\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"\"\n",
    "\n",
    "# Notebook at full width in the browser\n",
    "from IPython.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))\n",
    "\n",
    "import sys\n",
    "from urllib.request import urlretrieve\n",
    "from pathlib import Path\n",
    "from collections import defaultdict\n",
    "from abc import ABC, abstractmethod\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "matplotlib.rcParams[\"image.interpolation\"] = \"none\"\n",
    "matplotlib.rcParams['figure.figsize'] = (14, 10)\n",
    "import numpy as np\n",
    "from tifffile import imread, imwrite\n",
    "from tqdm.auto import tqdm\n",
    "import skimage\n",
    "import pandas as pd\n",
    "import scipy\n",
    "\n",
    "from stardist import fill_label_holes, random_label_cmap\n",
    "from stardist.plot import render_label\n",
    "from stardist.models import StarDist2D\n",
    "from stardist import _draw_polygons\n",
    "from csbdeep.utils import normalize\n",
    "\n",
    "import napari\n",
    "\n",
    "lbl_cmap = random_label_cmap()\n",
    "# Pretty tqdm progress bars \n",
    "! jupyter nbextension enable --py widgetsnbextension"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e690d63-437a-4850-9071-774d502f2ae3",
   "metadata": {},
   "source": [
    "Some utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a0c88d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_img_label(img, lbl, img_title=\"image\", lbl_title=\"label\", **kwargs):\n",
    "    fig, (ai,al) = plt.subplots(1,2, gridspec_kw=dict(width_ratios=(1,1)))\n",
    "    im = ai.imshow(img, cmap='gray', clim=(0,1))\n",
    "    ai.set_title(img_title)\n",
    "    ai.axis(\"off\")\n",
    "    al.imshow(render_label(lbl, img=.3*img, normalize_img=False, cmap=lbl_cmap))\n",
    "    al.set_title(lbl_title)\n",
    "    al.axis(\"off\")\n",
    "    plt.tight_layout()\n",
    "    \n",
    "def preprocess(X, Y, axis_norm=(0,1)):\n",
    "    # normalize channels independently\n",
    "    X = np.stack([normalize(x, 1, 99.8, axis=axis_norm) for x in tqdm(X, leave=True, desc=\"Normalize images\")])\n",
    "    # fill holes in labels\n",
    "    Y = np.stack([fill_label_holes(y) for y in tqdm(Y, leave=True, desc=\"Fill holes in labels\")])\n",
    "    return X, Y"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2a3c35e-328e-4557-b28c-770206268690",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## Inspect the dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41d42b15-b6a5-4663-b25a-0a272b205d43",
   "metadata": {},
   "source": [
    "For this exercise we will be working with a fluorenscence microscopy time-lapse of breast cancer cells with stained nuclei (SiR-DNA), originally from https://zenodo.org/record/4034976#.YwZRCJPP1qt.\n",
    "\n",
    "We will use a slightly modified version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71944856-44f7-454e-ad24-18e8cd0f6248",
   "metadata": {},
   "outputs": [],
   "source": [
    "base_path = Path(\"data/exercise1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72806388-6422-4e6c-8702-a7f5ad45a460",
   "metadata": {},
   "source": [
    "Load the dataset (images and tracking annotations) from disk into this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36b598b0-c764-4431-b2da-7cfa8e98a3ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.stack([imread(xi) for xi in sorted((base_path / \"images\").glob(\"*.tif\"))])\n",
    "y = np.stack([imread(yi) for yi in sorted((base_path / \"gt_tracking\").glob(\"*.tif\"))])\n",
    "assert len(x) == len(y)\n",
    "print(f\"Number of images: {len(x)}\")\n",
    "print(f\"Image shape: {x[0].shape}\")\n",
    "links = np.loadtxt(base_path / \"gt_tracking\" / \"man_track.txt\", dtype=int)\n",
    "links = pd.DataFrame(data=links, columns=[\"track_id\", \"from\", \"to\", \"parent_id\"])\n",
    "print(\"Links\")\n",
    "links[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aed16327-6787-4e2b-99f0-bdbb7e3f4c44",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = preprocess(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b984eb94-6c2b-4ff9-be3d-73c552a359fc",
   "metadata": {},
   "source": [
    "Visualize some images (by changing `idx`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08f65b4d-d3b0-41f6-8a75-be71fffd97df",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = 0\n",
    "plot_img_label(x[idx], y[idx])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c29f78c-11aa-4536-bd7f-814741d8140c",
   "metadata": {},
   "source": [
    "This is ok to take a glimpse, but a dynamic viewer would be much better. Let's use [napari](https://napari.org/tutorials/fundamentals/getting_started.html) for this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3630281e-96fa-4623-b9c9-674dc24a16fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "viewer = napari.Viewer()\n",
    "viewer.add_image(x, name=\"image\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b35e8d4-d7ca-4c61-b762-7546b76fbd26",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\"><h3>Napari in a jupyter notebook:</h3>\n",
    "    \n",
    "- To have napari working in a jupyter notebook, you need to use up-to-date versions of napari, pyqt and pyqt5, as is the case in the conda environments provided together with this exercise.\n",
    "- When you are coding and debugging, close the napari viewer with `viewer.close()` to avoid problems with the two event loops of napari and jupyter.\n",
    "- **If a cell is not executed (empty square brackets on the left of a cell) despite you running it, running it a second time right after will usually work.**\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "209d41af-1c6a-43c1-a10b-2b42501ff5ed",
   "metadata": {},
   "source": [
    "Let's add the ground truth annotations. Now we can easily explore how the cells move over time.\n",
    "\n",
    "If you zoom in, you will note that the dense annotations are not perfect segmentations, but rather circles placed roughly in the center of each nucleus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b305f8e8-a10b-4846-bd66-5c613f8c7b39",
   "metadata": {},
   "outputs": [],
   "source": [
    "def visualize_tracks(viewer, y, links=None, name=\"\"):\n",
    "    \"\"\"Utility function to visualize segmentation and tracks\"\"\"\n",
    "    max_label = max(links.max(), y.max()) if links is not None else y.max()\n",
    "    colorperm = np.random.default_rng(42).permutation((np.arange(1, max_label + 2)))\n",
    "    tracks = []\n",
    "    for t, frame in enumerate(y):\n",
    "        centers = skimage.measure.regionprops(frame)\n",
    "        for c in centers:\n",
    "            tracks.append([colorperm[c.label], t, int(c.centroid[0]), int(c.centroid[1])])\n",
    "    tracks = np.array(tracks)\n",
    "    tracks = tracks[tracks[:, 0].argsort()]\n",
    "    \n",
    "    graph = {}\n",
    "    if links is not None:\n",
    "        divisions = links[links[:,3] != 0]\n",
    "        for d in divisions:\n",
    "            if colorperm[d[0]] not in tracks[:, 0] or colorperm[d[3]] not in tracks[:, 0]:\n",
    "                continue\n",
    "            graph[colorperm[d[0]]] = [colorperm[d[3]]]\n",
    "\n",
    "    viewer.add_labels(y, name=f\"{name}_detections\")\n",
    "    viewer.layers[f\"{name}_detections\"].contour = 3\n",
    "    viewer.add_tracks(tracks, name=f\"{name}_tracks\", graph=graph)\n",
    "    return tracks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36755434-ed79-41c3-8d18-4658d0ac88bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "visualize_tracks(viewer, y, links.to_numpy(), \"ground_truth\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d6df6e2-332f-497c-a007-f1e126ffe325",
   "metadata": {
    "incorrectly_encoded_metadata": "jp-MarkdownHeadingCollapsed=true",
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## Exercise 1.1\n",
    "<div class=\"alert alert-block alert-info\"><h3>Exercise 1.1: Highlight the cell divisions</h3>\n",
    "\n",
    "The visualization of the ground truth tracks are useful to grasp this video, but it is still hard see the cell divisions. Given the dense annotations `y` and the track links `links`, write a function to create a new layer that highlights the pairs of daughter cells just after mitosis.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30861a2c-621e-4d03-a02f-c7a1b960e415",
   "metadata": {},
   "source": [
    "Expected outcome:<br>\n",
    "<figure style=\"display:inline-block\">\n",
    "    <img src=\"figures/prediv.png\" width=\"400\" />\n",
    "    <figcaption>frame t</figcaption>\n",
    "</figure>\n",
    "<figure style=\"display:inline-block\">\n",
    "    <img src=\"figures/postdiv.png\" width=\"400\" />\n",
    "    <figcaption>frame t+1</figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2a6b3ad-311f-429e-aab2-4530fb3b85f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def visualize_divisions(viewer, y, links):\n",
    "    \"\"\"Utility function to visualize divisions\"\"\"\n",
    "    \n",
    "    ### YOUR CODE HERE ###\n",
    "    \n",
    "    divisions = np.zeros_like(y)\n",
    "    viewer.add_labels(divisions, name=\"divisions\")\n",
    "    return divisions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db691d6e-b5ae-44c5-bd67-b7b8f18408f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "visualize_divisions(viewer, y, links.to_numpy());"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "835ca5c3-96b9-4229-81b4-37669bc3ba11",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## Object detection using a pre-trained neural network"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "763f0ec8-5b09-4369-9eb7-686234008150",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Load a pretrained stardist model, detect nuclei in one image and visualize them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8f156b6-2add-4974-ba6f-a4f813b26d50",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "idx = 0\n",
    "model = StarDist2D.from_pretrained(\"2D_versatile_fluo\")\n",
    "(detections, details), (prob, _) = model.predict_instances(x[idx], scale=(1, 1), return_predict=True)\n",
    "plot_img_label(x[idx], detections, lbl_title=\"detections\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85bff23f-7029-4594-ae2c-ead38652713d",
   "metadata": {},
   "source": [
    "Here we visualize in detail the polygons we have detected with StarDist. TODO some description on how StarDist works.\n",
    "\n",
    "<!-- Notice that each object comes with a center point, which we can use to compute pairwise euclidian distances between objects. -->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb767441-bab9-426f-83e7-5c6b7d094ad3",
   "metadata": {},
   "outputs": [],
   "source": [
    "coord, points, polygon_prob = details['coord'], details['points'], details['prob']\n",
    "plt.figure(figsize=(24,12))\n",
    "plt.subplot(121)\n",
    "plt.title(\"Predicted Polygons\")\n",
    "_draw_polygons(coord, points, polygon_prob, show_dist=True)\n",
    "plt.imshow(x[idx], cmap='gray'); plt.axis('off')\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.title(\"Object center probability\")\n",
    "plt.imshow(prob, cmap='magma'); plt.axis('off')\n",
    "plt.tight_layout()\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b06258ea-2cbe-43c6-84b9-e6e45e08082d",
   "metadata": {
    "incorrectly_encoded_metadata": "tags=[] jp-MarkdownHeadingCollapsed=true tags=[] jp-MarkdownHeadingCollapsed=true",
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## Exercise 1.2\n",
    "<div class=\"alert alert-block alert-info\"><h3>Exercise 1.2: Explore the parameters of cell detection</h3>\n",
    "\n",
    "Explore the following aspects of the detection algorithm:     \n",
    "- The `scale` parameter of the function `predict_instances` downscales (< 1) or upscales (> 1) the images by the given factor before feeding them to the neural network. How do the detections change if you adjust it?\n",
    "- Inspect false positive and false negative detections. Do you observe patterns?\n",
    "- So far we have used a StarDist off the shelf. Luckily, we also have a StarDist model that was trained on a subset of this breast cancer cells dataset (from https://zenodo.org/record/4034976#.Yv-aNPFBzao). Load it with `model = StarDist2D(None, name=\"stardist_breast_cancer\", basedir=\"models\")` and qualitatively observe differences.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8faf20a8-b085-4e1f-b3ee-f0af8798d4c2",
   "metadata": {},
   "source": [
    "Detect centers and segment nuclei in all images of the time lapse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bc2e533-58b9-4b5b-9c61-adb9cdd8a245",
   "metadata": {},
   "outputs": [],
   "source": [
    "scale = (1.0, 1.0)\n",
    "pred = [model.predict_instances(xi, show_tile_progress=False, scale=scale)\n",
    "              for xi in tqdm(x)]\n",
    "detections = [xi[0] for xi in pred]\n",
    "detections = np.stack([skimage.segmentation.relabel_sequential(d)[0] for d in detections])  # ensure that label ids are contiguous and start at 1 for each frame \n",
    "centers = [xi[1][\"points\"] for xi in pred]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf114a1f-2d24-4976-8cdd-953d93c31469",
   "metadata": {},
   "source": [
    "Visualize the dense detections. Note that they are still not linked and therefore randomly colored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f997dd09-3121-402b-b7b9-9e5498bad2a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "viewer = napari.viewer.current_viewer()\n",
    "if viewer:\n",
    "    viewer.close()\n",
    "viewer = napari.Viewer()\n",
    "viewer.add_image(x)\n",
    "viewer.add_labels(detections, name=f\"detections_scale_{scale}\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2aed356a-d8e7-4b0d-a89c-cc8b61269f0c",
   "metadata": {},
   "source": [
    "We see that the number of detections increases over time, corresponding to the cells that insert the field of view from below during the video."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff23798d-57d7-4a18-a00c-9b7d5982e235",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,6))\n",
    "plt.bar(range(len(centers)), [len(xi) for xi in centers])\n",
    "plt.title(f\"Number of detections in each frame (scale={scale})\")\n",
    "plt.xticks(range(len(centers)))\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c033125-d4f5-48bb-9b2c-47aa7c8e5ca7",
   "metadata": {
    "incorrectly_encoded_metadata": "tags=[] jp-MarkdownHeadingCollapsed=true"
   },
   "source": [
    "## Checkpoint 1\n",
    "<div class=\"alert alert-block alert-success\"><h3>Checkpoint 1: We have good detections, now on to the linking.</h3></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b9f5a9a-dfe4-4568-a85d-a1dfc3eba19a",
   "metadata": {
    "incorrectly_encoded_metadata": "tags=[] jp-MarkdownHeadingCollapsed=true",
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## Greedy linking by nearest neighbor\n",
    "\n",
    "TODO write introduction text."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbc3769f-6135-4b42-a779-610353fa7cc1",
   "metadata": {},
   "source": [
    "## Exercise 1.3\n",
    "<div class=\"alert alert-block alert-info\"><h3>Exercise 1.3: Write a function that computes pairwise euclidian distances given two lists of points.</h3></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d1ff258-73ba-4fe8-8636-c6f7e65fbfb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def euclidian_distance(points0, points1):\n",
    "    dists = np.zeros((len(points0), len(points1)))\n",
    "    \n",
    "    ### YOUR CODE HERE ###\n",
    "    \n",
    "    return dists"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8141410c-da23-4fbc-8755-b14b0f79e8e4",
   "metadata": {
    "tags": []
   },
   "source": [
    "Here are two (almost random ;)) lists of points to test your function on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ec05e1f-436d-4309-9cd4-1f47a30ce189",
   "metadata": {},
   "outputs": [],
   "source": [
    "green_points = np.load(\"points.npz\")[\"green\"]\n",
    "cyan_points = np.load(\"points.npz\")[\"cyan\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec69da1e-cab4-473c-85c1-c45fb1e2ba16",
   "metadata": {},
   "outputs": [],
   "source": [
    "%time dists = pairwise_euclidian_distance(green_points, cyan_points)\n",
    "assert np.allclose(dists, np.load(\"points.npz\")[\"dists_green_cyan\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c8a72df-3715-4a24-85a9-1f87a86e85ef",
   "metadata": {},
   "source": [
    "## Exercise 1.4\n",
    "<div class=\"alert alert-block alert-info\"><h3>Exercise 1.4: Write a function that greedily extracts a nearest neighbors assignment given a cost matrix.</h3></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6880a24a-7942-4cb8-84fa-ddab9396c794",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nearest_neighbor(cost_matrix):\n",
    "    \"\"\"Greedy nearest neighbor assignment.\n",
    "    \n",
    "    Each point in both sets can only be assigned once. \n",
    "    \n",
    "    Args:\n",
    "\n",
    "        cost_matrix: m x n matrix with pairwise linking costs of two sets of points.\n",
    "\n",
    "    Returns:\n",
    "\n",
    "        Determined matches as tuple of lists (ids_of_rows, ids_of_columns).\n",
    "    \"\"\"\n",
    "\n",
    "    ids_from = []\n",
    "    ids_to = []\n",
    "    \n",
    "    ### YOUR CODE HERE ###\n",
    "    \n",
    "    return np.array(ids_from), np.array(ids_to)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "883a71e2-2e56-4309-8831-09aeb0f0bd5c",
   "metadata": {},
   "source": [
    "Test your implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "888d1056-6502-4f55-9c6c-237308faf14b",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_matrix = np.array([\n",
    "    [9, 2, 9],\n",
    "    [9, 9, 9],\n",
    "    [1, 9, 9],\n",
    "    [9, 3, 9],\n",
    "])\n",
    "idx_from, idx_to = nearest_neighbor(test_matrix)\n",
    "assert np.all(idx_from == [2, 0, 1])\n",
    "assert np.all(idx_to == [0, 1, 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bca8c7dc-e44e-4b5c-b534-b4a4ef3130fe",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Exercise 1.5\n",
    "<div class=\"alert alert-block alert-info\"><h3>Exercise 1.5: Complete a thresholded nearest neighbor linker using your functions from exercises 1.3 and 1.4.</h3>\n",
    "\n",
    "You have to write two methods:\n",
    "    \n",
    "- Method 1 (`linking_cost_function`): Given dense detections in a pair of frames, extract their centroids and calculate pairwise euclidian distances between them. \n",
    "- Method 2 (`_link_two_frames`): For each detection in frame $t$, find the nearest neighbor in frame $t+1$ given the cost matrix. If the distance is below a threshold $\\tau$, link the two objects. Explore different values of threshold $\\tau$.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd5d1a2d-5ab1-4f7d-8408-64b54ebfa7e7",
   "metadata": {},
   "source": [
    "Here you are seeing an abstract base class (`ABC`) for linking detections in a video with some local frame-by-frame algorithm.\n",
    "\n",
    "The class already comes with some useful methods that you won't have to worry about, such as iterating over frames, visualizing linked results as well as sanity checks of inputs.\n",
    "\n",
    "There are two abstract methods (\"gaps\") in `FrameByFrameLinker`:\n",
    "- `linking_cost_function`\n",
    "- `_link_two_frames`\n",
    "\n",
    "In the exercises 1.5 - 1.8, you will make different subclasses of `FrameByFrameLinker`, in which it will be your job to write these two methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d96c279f-5e0f-47f5-b35f-4d275a486161",
   "metadata": {},
   "outputs": [],
   "source": [
    "class FrameByFrameLinker(ABC):\n",
    "    \"\"\"Abstract base class for linking detections by considering pairs of adjacent frames.\"\"\"\n",
    "    \n",
    "    def link(self, detections, images=None):\n",
    "        \"\"\"Links detections in t frames.\n",
    "        \n",
    "        Args:\n",
    "        \n",
    "            detections:\n",
    "            \n",
    "                List of t numpy arrays of shape (x,y) with contiguous label ids. Background = 0.\n",
    "                \n",
    "            images (optional):\n",
    "            \n",
    "                List of t numpy arrays of shape (x,y).\n",
    "        \n",
    "        Returns:\n",
    "        \n",
    "            List of t linking dictionaries, each containing:\n",
    "                \"links\": Tuple of lists (ids frame t, ids frame t+1),\n",
    "                \"births\": List of ids,\n",
    "                \"deaths\": List of ids.\n",
    "            Ids are one-based, 0 is reserved for background.\n",
    "        \"\"\"\n",
    "        if images is not None:\n",
    "            assert len(images) == len(detections)\n",
    "        else:\n",
    "            images = [None] * len(detections)\n",
    "\n",
    "        links = []\n",
    "        for i in tqdm(range(len(images) - 1), desc=\"Linking\"):\n",
    "            detections0 = detections[i]\n",
    "            detections1 = detections[i+1]\n",
    "            self._assert_relabeled(detections0)\n",
    "            self._assert_relabeled(detections1)\n",
    "            \n",
    "            cost_matrix = self.linking_cost_function(detections0, detections1, images[i], images[i+1])\n",
    "            li = self._link_two_frames(cost_matrix)\n",
    "            self._assert_links(links=li, time=i, detections0=detections0, detections1=detections1) \n",
    "            links.append(li)\n",
    "            \n",
    "        return links\n",
    "\n",
    "    @abstractmethod\n",
    "    def linking_cost_function(self, detections0, detections1, image0=None, image1=None):\n",
    "        \"\"\"Calculate features for each detection and extract pairwise costs.\n",
    "        \n",
    "        To be overwritten in subclass.\n",
    "        \n",
    "        Args:\n",
    "        \n",
    "            detections0: image with background 0 and detections 1, ..., m\n",
    "            detections1: image with backgruond 0 and detections 1, ..., n\n",
    "            image0 (optional): image corresponding to detections0\n",
    "            image1 (optional): image corresponding to detections1\n",
    "            \n",
    "        Returns:\n",
    "        \n",
    "            m x n cost matrix \n",
    "        \"\"\"\n",
    "        pass\n",
    "    \n",
    "    @abstractmethod\n",
    "    def _link_two_frames(self, cost_matrix):\n",
    "        \"\"\"Link two frames.\n",
    "        \n",
    "        To be overwritten in subclass.\n",
    "\n",
    "        Args:\n",
    "\n",
    "            cost_matrix: m x n matrix\n",
    "\n",
    "        Returns:\n",
    "        \n",
    "            Linking dictionary:\n",
    "                \"links\": Tuple of lists (ids frame t, ids frame t+1),\n",
    "                \"births\": List of ids,\n",
    "                \"deaths\": List of ids.\n",
    "            Ids are one-based, 0 is reserved for background.\n",
    "        \"\"\"\n",
    "        pass\n",
    "\n",
    "    def relabel_detections(self, detections, links):\n",
    "        \"\"\"Relabel dense detections according to computed links, births and deaths.\n",
    "        \n",
    "        Args:\n",
    "        \n",
    "            detections: \n",
    "                 \n",
    "                 List of t numpy arrays of shape (x,y) with contiguous label ids. Background = 0.\n",
    "                 \n",
    "            links:\n",
    "                \n",
    "                List of t linking dictionaries, each containing:\n",
    "                    \"links\": Tuple of lists (ids frame t, ids frame t+1),\n",
    "                    \"births\": List of ids,\n",
    "                    \"deaths\": List of ids.\n",
    "                Ids are one-based, 0 is reserved for background.\n",
    "        \"\"\"\n",
    "        detections = detections.copy()\n",
    "        \n",
    "        assert len(detections) - 1 == len(links)\n",
    "        self._assert_relabeled(detections[0])\n",
    "        out = [detections[0]]\n",
    "        n_tracks = out[0].max()\n",
    "        lookup_tables = [{i: i for i in range(1, out[0].max() + 1)}]\n",
    "\n",
    "        for i in tqdm(range(len(links)), desc=\"Recoloring detections\"):\n",
    "            (ids_from, ids_to) = links[i][\"links\"]\n",
    "            births = links[i][\"births\"]\n",
    "            deaths = links[i+1][\"deaths\"] if i+1 < len(links) else []\n",
    "            new_frame = np.zeros_like(detections[i+1])\n",
    "            self._assert_relabeled(detections[i+1])\n",
    "            \n",
    "            lut = {}\n",
    "            for _from, _to in zip(ids_from, ids_to):\n",
    "                # Copy over ID\n",
    "                new_frame[detections[i+1] == _to] = lookup_tables[i][_from]\n",
    "                lut[_to] = lookup_tables[i][_from]\n",
    "\n",
    "            \n",
    "            # Start new track for birth tracks\n",
    "            for b in births:\n",
    "                if b in deaths:\n",
    "                    continue\n",
    "                \n",
    "                n_tracks += 1\n",
    "                lut[b] = n_tracks\n",
    "                new_frame[detections[i+1] == b] = n_tracks\n",
    "                \n",
    "            # print(lut)\n",
    "            lookup_tables.append(lut)\n",
    "            out.append(new_frame)\n",
    "                \n",
    "        return np.stack(out)\n",
    "\n",
    "    def _assert_links(self, links, time, detections0, detections1):\n",
    "        if len(links[\"links\"][0]) != len(links[\"links\"][1]):\n",
    "            raise RuntimeError(\"Format of links['links'] not correct.\")\n",
    "            \n",
    "        if sorted([*links[\"links\"][0], *links[\"deaths\"]]) != list(range(1, len(np.unique(detections0)))):\n",
    "            raise RuntimeError(f\"Some detections in frame {time} are not properly assigned as either linked or death.\")\n",
    "            \n",
    "        if sorted([*links[\"links\"][1], *links[\"births\"]]) != list(range(1, len(np.unique(detections1)))):\n",
    "            raise RuntimeError(f\"Some detections in frame {time + 1} are not properly assigned as either linked or birth.\")\n",
    "            \n",
    "        for b in links[\"births\"]:\n",
    "            if b in links[\"links\"][1]:\n",
    "                raise RuntimeError(f\"Links frame {time+1}: Detection {b} marked as birth, but also linked.\")\n",
    "        \n",
    "        for d in links[\"deaths\"]:\n",
    "            if d in links[\"links\"][0]:\n",
    "                raise RuntimeError(f\"Links frame {time}: Detection {d} marked as death, but also linked.\")\n",
    "        \n",
    "        \n",
    "    def _assert_relabeled(self, x):\n",
    "        if x.min() < 0:\n",
    "            raise ValueError(\"Negative ID in detections.\")\n",
    "        if x.min() == 0:\n",
    "            n = x.max() + 1\n",
    "        else:\n",
    "            n = x.max()\n",
    "        if n != len(np.unique(x)):\n",
    "            raise ValueError(\"Detection IDs are not contiguous.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c72b377-bbcc-4a1f-991b-900b6f630a11",
   "metadata": {},
   "source": [
    "Hints:\n",
    "- Check out `skimage.measure.regionprops`.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89deee7b-881c-467d-ab79-78f29b207fa3",
   "metadata": {},
   "outputs": [],
   "source": [
    "class NearestNeighborLinkerEuclidian(FrameByFrameLinker):\n",
    "    \"\"\".\n",
    "    \n",
    "    Args:\n",
    "    \n",
    "        threshold (float): Maximum euclidian distance for linking.\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, threshold=sys.float_info.max, *args, **kwargs):\n",
    "        self.threshold = threshold\n",
    "        super().__init__(*args, **kwargs)\n",
    "    \n",
    "    def linking_cost_function(self, detections0, detections1, image0=None, image1=None):\n",
    "        \"\"\" Get centroids from detections and compute pairwise euclidian distances.\n",
    "                \n",
    "        Args:\n",
    "        \n",
    "            detections0: image with background 0 and detections 1, ..., m\n",
    "            detections1: image with backgruond 0 and detections 1, ..., n\n",
    "            \n",
    "        Returns:\n",
    "        \n",
    "            m x n cost matrix \n",
    "        \"\"\"\n",
    "        ### YOUR CODE HERE ###\n",
    "        dists = np.zeros((detections0.max(), detections1.max()))\n",
    "        \n",
    "        return dists\n",
    "    \n",
    "    def _link_two_frames(self, cost_matrix):\n",
    "        \"\"\"Greedy nearest neighbor assignment.\n",
    "\n",
    "        Each point in both sets can only be assigned once. \n",
    "\n",
    "        Args:\n",
    "\n",
    "            cost_matrix: m x n matrix containing pairwise linking costs of two sets of points.\n",
    "\n",
    "        Returns:\n",
    "            Linking dictionary:\n",
    "                \"links\": Tuple of lists (ids frame t, ids frame t+1),\n",
    "                \"births\": List of ids,\n",
    "                \"deaths\": List of ids.\n",
    "            Ids are one-based, 0 is reserved for background.\n",
    "        \"\"\"\n",
    "        \n",
    "        min_objs = min(cost_matrix.shape[0], cost_matrix.shape[1])\n",
    "        ids_from = np.arange(min_objs)\n",
    "        ids_to = np.arange(min_objs)\n",
    "        births = np.arange(min_objs, cost_matrix.shape[1])\n",
    "        deaths = np.arange(min_objs, cost_matrix.shape[0])\n",
    "        \n",
    "        ### YOUR CODE HERE (REPLACE THE DUMMY INITIALIZATIONS ABOVE) ###\n",
    "        \n",
    "                            \n",
    "        # Account for +1 offset of the dense labels\n",
    "        ids_from += 1\n",
    "        ids_to += 1\n",
    "        births += 1\n",
    "        deaths += 1\n",
    "        \n",
    "        links = {\"links\": (ids_from, ids_to), \"births\": births, \"deaths\": deaths}\n",
    "        return links"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d005234b-9c4b-4c60-8656-00d0a92790fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# nn_linker = NearestNeighborLinkerEuclidian(threshold=1000) # Explore different values of `threshold`\n",
    "nn_linker = NearestNeighborLinkerEuclidian(threshold=50) # Solution param\n",
    "nn_links = nn_linker.link(detections)\n",
    "nn_tracks = nn_linker.relabel_detections(detections, nn_links)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40291085-f94c-4925-857b-00c5162829b7",
   "metadata": {},
   "source": [
    "Visualize results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76932eb8-0186-4777-8818-8e976940f965",
   "metadata": {},
   "outputs": [],
   "source": [
    "viewer = napari.viewer.current_viewer()\n",
    "if viewer:\n",
    "    viewer.close()\n",
    "viewer = napari.Viewer()\n",
    "viewer.add_image(x)\n",
    "visualize_tracks(viewer, nn_tracks, name=\"nn\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5ab42b7-ef07-4476-bc7d-d329ea91636b",
   "metadata": {},
   "source": [
    "## Checkpoint 2\n",
    "<div class=\"alert alert-block alert-success\"><h3>Checkpoint 2: We built a basic tracking algorithm from scratch :).</h3></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f241c4fa-e205-4de4-b9c6-7c0237459f9d",
   "metadata": {
    "incorrectly_encoded_metadata": "tags=[] jp-MarkdownHeadingCollapsed=true"
   },
   "source": [
    "## Exercise 1.6\n",
    "<div class=\"alert alert-block alert-info\"><h3>Exercise 1.6: Estimate the global drift of the data</h3>\n",
    "\n",
    "We can observe that all cells move upwards with an approximately constant displacement in each timestep. Write a slightly modified version of `NearestNeighborLinkerEuclidian` with a slightly modified `linking_cost_function` that accounts for this.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16e7f0eb-19f1-42e8-baae-52ddb0a5c279",
   "metadata": {},
   "outputs": [],
   "source": [
    "class NearestNeighborLinkerDriftCorrection(NearestNeighborLinkerEuclidian):\n",
    "    \"\"\".\n",
    "    \n",
    "    Args:\n",
    "        \n",
    "        drift: tuple of (x,y) drift correction per frame.\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, drift, *args, **kwargs):\n",
    "        self.drift = np.array(drift)\n",
    "        super().__init__(*args, **kwargs)\n",
    "        \n",
    "    def linking_cost_function(self, detections0, detections1, image0=None, image1=None):\n",
    "        \"\"\" Get centroids from detections and compute pairwise euclidian distances with drift correction.\n",
    "                \n",
    "        Args:\n",
    "        \n",
    "            detections0: image with background 0 and detections 1, ..., m\n",
    "            detections1: image with backgruond 0 and detections 1, ..., n\n",
    "            \n",
    "        Returns:\n",
    "        \n",
    "            m x n cost matrix \n",
    "        \"\"\"\n",
    "        \n",
    "        ### YOUR CODE HERE\n",
    "        \n",
    "        dists = np.zeros((detections0.max(), detections1.max()))\n",
    "        return dists"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a706de6-da8f-4503-870b-86fd4aca123b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Explore different values of `threshold` and `drift`\n",
    "# drift_linker = NearestNeighborLinkerDriftCorrection(threshold=1000, drift=(0, 0))\n",
    "drift_linker = NearestNeighborLinkerDriftCorrection(threshold=50, drift=(-20, 0)) # SOLUTION params\n",
    "drift_links = drift_linker.link(detections)\n",
    "drift_tracks = drift_linker.relabel_detections(detections, drift_links)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f36e7a0f-fe86-477f-a282-3b296e35c302",
   "metadata": {},
   "source": [
    "Visualize results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbf2c9e3-6dbc-4f4a-b66b-44f7852396ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "viewer = napari.viewer.current_viewer()\n",
    "if viewer:\n",
    "    viewer.close()\n",
    "viewer = napari.Viewer()\n",
    "viewer.add_image(x)\n",
    "visualize_tracks(viewer, drift_tracks, name=\"drift\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c96e9e0-0fcc-4f2b-b635-3797df6e8ff9",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Optimal frame-by-frame matching (*Linear assignment problem* or *Weighted bipartite matching*)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c73cb9b-81c2-4d7b-8a53-3ab2497cb995",
   "metadata": {
    "incorrectly_encoded_metadata": "jp-MarkdownHeadingCollapsed=true",
    "tags": []
   },
   "source": [
    "## Exercise 1.7\n",
    "<div class=\"alert alert-block alert-info\"><h3>Exercise 1.7: Perform optimal frame-by-frame linking</h3>\n",
    "\n",
    "Set up the cost matrix such that you can use [`scipy.optimize.linear_sum_assignment`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linear_sum_assignment.html) to solve the matching problem in the bipartite graph.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1bd64d6-7428-41dc-b483-47bad2ae4627",
   "metadata": {},
   "source": [
    "TODO write intro text.\n",
    "\n",
    "\n",
    "TODO insert image for bipartite matching\n",
    "\n",
    "<img src=\"figures/LAP_cost_matrix.png\" width=\"300\"/>\n",
    "\n",
    "\n",
    "from Jaqaman, Khuloud, et al. \"Robust single-particle tracking in live-cell time-lapse sequences.\" Nature Methods (2008)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f8db520-59f2-4d81-ab57-9f545fba518e",
   "metadata": {},
   "outputs": [],
   "source": [
    "class BipartiteMatchingLinker(FrameByFrameLinker):\n",
    "    \"\"\".\n",
    "    \n",
    "    Args:\n",
    "        threshold (float): Maximum euclidian distance for linking.\n",
    "        drift: tuple of (x,y) drift correction per frame.\n",
    "        birth_cost_factor (float): Multiply factor with maximum entry in cost matrix.\n",
    "        death_cost_factor (float): Multiply factor with maximum entry in cost matrix.\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(\n",
    "        self,\n",
    "        threshold=sys.float_info.max,\n",
    "        drift=(0,0),\n",
    "        birth_cost_factor=1.05,\n",
    "        death_cost_factor=1.05,\n",
    "        *args,\n",
    "        **kwargs\n",
    "    ):\n",
    "        self.threshold = threshold\n",
    "        self.drift = np.array(drift)\n",
    "        self.birth_cost_factor = birth_cost_factor\n",
    "        self.death_cost_factor = death_cost_factor\n",
    "        \n",
    "        super().__init__(*args, **kwargs)\n",
    "        \n",
    "    def linking_cost_function(self, detections0, detections1, image0=None, image1=None):\n",
    "        \"\"\" Get centroids from detections and compute pairwise euclidian distances with drift correction.\n",
    "                \n",
    "        Args:\n",
    "        \n",
    "            detections0: image with background 0 and detections 1, ..., m\n",
    "            detections1: image with backgruond 0 and detections 1, ..., n\n",
    "            \n",
    "        Returns:\n",
    "        \n",
    "            m x n cost matrix \n",
    "        \"\"\"\n",
    "        dists = np.zeros((detections0.max(), detections1.max()))\n",
    "        return dists\n",
    "    \n",
    "    def _link_two_frames(self, cost_matrix):\n",
    "        \"\"\"Weighted bipartite matching with square matrix from Jaqaman et al (2008).\n",
    "\n",
    "        Args:\n",
    "\n",
    "            cost_matrix: m x n matrix.\n",
    "\n",
    "        Returns:\n",
    "        \n",
    "            Linking dictionary:\n",
    "                \"links\": Tuple of lists (ids frame t, ids frame t+1),\n",
    "                \"births\": List of ids,\n",
    "                \"deaths\": List of ids.\n",
    "            Ids are one-based, 0 is reserved for background.\n",
    "        \"\"\"\n",
    "        \n",
    "        cost_matrix = cost_matrix.copy().astype(float)\n",
    "        b = self.birth_cost_factor * min(self.threshold, cost_matrix.max())\n",
    "        d = self.death_cost_factor * min(self.threshold, cost_matrix.max())\n",
    "        no_link = max(cost_matrix.max(), max(b, d)) * 1e9\n",
    "        \n",
    "        \n",
    "        min_objs = min(cost_matrix.shape[0], cost_matrix.shape[1])\n",
    "        ids_from = np.arange(min_objs)\n",
    "        ids_to = np.arange(min_objs)\n",
    "        births = np.arange(min_objs, cost_matrix.shape[1])\n",
    "        deaths = np.arange(min_objs, cost_matrix.shape[0])\n",
    "        \n",
    "        ### YOUR CODE HERE (REPLACE THE DUMMY INITIALIZATIONS FOR THE RETURN VARIABLES ABOVE) ###\n",
    "                        \n",
    "        # Account for +1 offset of the dense labels\n",
    "        ids_from += 1\n",
    "        ids_to += 1\n",
    "        births += 1\n",
    "        deaths += 1\n",
    "        \n",
    "        links = {\"links\": (ids_from, ids_to), \"births\": births, \"deaths\": deaths}\n",
    "        return links"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "697622ea-e7fd-4322-b538-a754106f2b0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "bm_linker = BipartiteMatchingLinker(threshold=50, drift=(-20, 0), birth_cost_factor=1.05, death_cost_factor=1.05)\n",
    "bm_links = bm_linker.link(detections)\n",
    "bm_tracks = bm_linker.relabel_detections(detections, bm_links)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e92ec76-0bda-4565-9c3b-34ff43c99ea6",
   "metadata": {},
   "outputs": [],
   "source": [
    "viewer = napari.viewer.current_viewer()\n",
    "if viewer:\n",
    "    viewer.close()\n",
    "viewer = napari.Viewer()\n",
    "viewer.add_image(x)\n",
    "visualize_tracks(viewer, bm_tracks, name=\"bm\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d27e1a02-739c-43de-ae50-c893d9e1a71e",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Other suitable features for linking cost function"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9385bdf3-8158-4d6e-a4c7-1f29e32df976",
   "metadata": {},
   "source": [
    "## Exercise 1.8\n",
    "\n",
    "<div class=\"alert alert-block alert-info\"><h3>Exercise 1.8: Explore different features for assigment problem</h3>\n",
    "\n",
    "Explore solving the assignment problem based different features and cost functions.\n",
    "For example:\n",
    "- Different morphological properties of detections (e.g. using `skimage.measure.regionprops`).\n",
    "- Extract texture features from the images, e.g. mean intensity for each detection.\n",
    "- Pairwise *Intersection over Union (IoU)* of detections.\n",
    "- ...\n",
    "\n",
    "Feel free to share features that improved the results with the class :).    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a57ae83-0833-41a8-91dd-fe676f82aa24",
   "metadata": {},
   "outputs": [],
   "source": [
    "class YourLinker(BipartiteMatchingLinker):\n",
    "    \n",
    "    def __init__(self, *args, **kwargs):\n",
    "        super().__init__(*args, **kwargs)\n",
    "    \n",
    "    def linking_cost_function(self, detections0, detections1, image0=None, image1=None):\n",
    "        \"\"\" Your very smart cost function for frame-by-frame linking.\n",
    "                \n",
    "        Args:\n",
    "        \n",
    "            detections0: image with background 0 and detections 1, ..., m\n",
    "            detections1: image with backgruond 0 and detections 1, ..., n\n",
    "            image0 (optional): image corresponding to detections0\n",
    "            image1 (optional): image corresponding to detections1\n",
    "            \n",
    "        Returns:\n",
    "        \n",
    "            m x n cost matrix \n",
    "        \"\"\"\n",
    "        return np.zeros((detections0.max(), detections1.max()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbd60b85-4bde-4e1a-bb8a-0a35c1514b5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "your_linker = YourLinker()\n",
    "your_links = your_linker.link(detections)\n",
    "your_tracks = your_linker.relabel_detections(detections, your_links)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa4a293a-4026-4bfb-b0c0-cfaeb5b2ca70",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,py:percent"
  },
  "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}