sessions/day04-tracking-and-time-series/02_tracking/Tracking_worksheet_STUDENT_VERSION.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Object tracking in Python\n",
    "**Stephen Cross (Wolfson Bioimaging Facility, University of Bristol, 2019).  Email: stephen.cross@bristol.ac.uk**\n",
    "\n",
    "In this session we will build up a workflow to track pre-detected objects through multiple frames of a timeseries.  \n",
    "\n",
    "The images we will be using are of fluorescently-labelled nuclei and are provided by the [Cell Tracking Challenge](http://celltrackingchallenge.net/).  For more information on the images, see the paper \"Recruitment of Oct4 Protein to UV-Damaged Chromatin in Embryonic Stem Cells\", Eva Bártová ,Gabriela Šustáčková,Lenka Stixová,Stanislav Kozubek,Soňa Legartová,Veronika Foltánková (2011) PLOS ONE.\n",
    "\n",
    "The session is broken down into eight exercises, each of which will introduce a new step to the workflow.  Typically an exercise will require creation of a new Python function to achieve a specific task (e.g. calculating the cost of linking two points).  The format of each function is laid out in the comments above that function.  It'll be up to you to replace the comment \"### CODE GOES HERE ###\" with the relevant code.\n",
    "\n",
    "At the end of each exercise is a pre-prepared block of code (\"### TEST SPACE ###\") that can be used to test the function is working correctly.  The code in this box doesn't need to be altered to work."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 0 - Getting started\n",
    "Before we start, we need to make sure this Jupyter Notebook has access to all the functions we're going to use later on.  \n",
    "\n",
    "For the sake of keeping the output on the notebook neat, we'll also set Numpy and Pandas to only display numbers to 2 decimal places (they're still working with higher precision values).  \n",
    "\n",
    "Finally, we'll also create a variable to represent \"infinity\".  This will be used later on during the object linking steps.\n",
    "\n",
    "### Aim\n",
    "- Import any functions we will use later on\n",
    "- Set Numpy and Pandas to display numbers to 2 decimal places\n",
    "- Create an \"infinity\" variable, which we will use during tracking\n",
    "- If not already completed, download and extract the images (this will all be done automatically)\n",
    "\n",
    "### Notes\n",
    "- This exercise doesn't require addition of any code; all you need do is run the code blocks to check all the libraries import correctly. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Importing the libraries into this Notebook\n",
    "import math\n",
    "import sys\n",
    "import util\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from scipy.optimize import linear_sum_assignment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setting parameters for Jupyter output\n",
    "np.set_printoptions(precision=2,threshold=sys.maxsize)\n",
    "pd.options.display.float_format = \"{:.2f}\".format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating a global variable for \"infinity\".  This will be needed later on when calculating assignments.\n",
    "inff = 1000000000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Getting the data (courtesy of Chas Nelson!)\n",
    "import urllib.request\n",
    "import io\n",
    "import os\n",
    "import zipfile\n",
    "\n",
    "from tqdm import tqdm_notebook, tnrange\n",
    "\n",
    "url = 'http://data.celltrackingchallenge.net/training-datasets/Fluo-N2DH-GOWT1.zip'\n",
    "\n",
    "if os.path.isfile (\"./assets/images/Fluo-N2DH-GOWT1/01/t001.tif\"):\n",
    "    print(\"Images already downloaded\")\n",
    "else:\n",
    "    with urllib.request.urlopen(url) as response:\n",
    "        print(\"Downloading...\")\n",
    "        length = int(response.getheader('content-length'))\n",
    "        chunk = max(4096, length//9999)\n",
    "    \n",
    "        buffer = io.BytesIO()\n",
    "        size = 0\n",
    "        for b in tnrange(length//chunk + 1):\n",
    "            block = response.read(chunk)\n",
    "            if not block:\n",
    "                print(\"Finished reading after {0}% of file.\".format(size/length))\n",
    "            buffer.write(block)\n",
    "            size = size + len(block)\n",
    "        print(\"Finished reading file.\")\n",
    "    \n",
    "        print(\"Unzipping... \",end=\"\")\n",
    "        zf = zipfile.ZipFile(buffer)\n",
    "        os.makedirs('./assets/images/',exist_ok=True)\n",
    "        zf.extractall(path='./assets/images')\n",
    "        print(\"Complete\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1 - Loading coordinates and visualising\n",
    "Our key aim this afternoon is to learn about tracking, so we don't want to spend half our time detecting objects to be tracked.  As such, you'll find a pre-prepared CSV file with object coordinates at *../assets/UntrackedCoordinates.csv*.  This file has 6 columns: Point ID, x-centroid, y-centroid, timepoint, 2D area and track ID (currently set to 0).\n",
    "\n",
    "We could also sink a couple of hours into creating a script to draw our object coordinates and tracks, but let's not.  For now, I've created a separate script which will load these files for us - this is in the \"util.py\" file we just imported.  This script also has functions to display the spots on top of the images, so we can check how our tracking is doing.\n",
    "\n",
    "### Aims\n",
    "- Load object coordinates from *../assets/UntrackedCoordinates.csv* into a Pandas dataframe\n",
    "- Load the timeseries images corresponding to the coordinates into a 3D Numpy array\n",
    "- Display the first few lines of the coordinate dataframe\n",
    "- Visualise the coordinates as an overlay on the timeseries images\n",
    "\n",
    "### Notes\n",
    "- This exercise also doesn't require any code to be written.  It's more about checking the relevant example files can be loaded and the utility functions run as expected.\n",
    "- If you get an error from \"%matplotlib notebook\", try changing it to \"%matplotlib widget\" then restarting the kernel (you'll also need to rerun previously-executed cells)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%html\n",
    "<style>\n",
    ".output_wrapper button.btn.btn-default,\n",
    ".output_wrapper .ui-dialog-titlebar {\n",
    "  display: none;\n",
    "}\n",
    "</style>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "\n",
    "# Loading image stack\n",
    "path = \"./assets/images/Fluo-N2DH-GOWT1/01/\"\n",
    "images = util.load_images(path);\n",
    "\n",
    "# Loading coordinates\n",
    "path = \"./assets/UntrackedCoordinates.csv\"\n",
    "coords = util.load_coordinates(path);\n",
    "\n",
    "# Displaying the first 10 lines of the dataframe\n",
    "print(coords.shape)\n",
    "print(coords[450:480])\n",
    "print(\"\")\n",
    "\n",
    "# Adding track renders\n",
    "util.show_overlay(images,coords,False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2 - Extracting coordinates for a specific timepoint\n",
    "We now have the full list of coordinates, but for tracking we're often going to want to access just those from a specific frame.  The first function we'll create will take the full list of coordinates and return a list of rows corresponding to only those points in the specified frame.\n",
    "\n",
    "### Aims\n",
    "- Create a function to get a list of row indices corresponding to coordinates in a specific frame\n",
    "- Arguments: The coordinates dataframe and the frame number as arguments\n",
    "- Return: A list of row indices locating the relevant points in the coordinates dataframe\n",
    "\n",
    "### Notes\n",
    "- This is the first exercise where we need to add in some code.  This exercise also follows the standard format that all remaining exercises will use.  In the next box you will need to replace the comment \"### CODE GOES HERE ###\" with the code for that function.  Above the function you'll find a description of the function format (i.e. the inputs and outputs).  Once the code is ready, you can test it with the next box (\"### TEST SPACE ###\").  The test space should display the row indices for the first couple of frames.  \n",
    "- If you wish to further verify these indices are correct, you could open the coordinates csv file in a spreadsheet viewer (e.g. Excel).\n",
    "- You'll likely want to use Panda's \"index\" function.  This returns the indices for the rows matching the provided condition.  For example, to get the row indices for all coordinates in frames after frame 7, we could use the following:\n",
    "```Python\n",
    "    coord_rows_over_f7 = coords.index[coords.FRAME > 7]\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to get the row indices for coordinates present in a specific frame.\n",
    "#\n",
    "# Args:\n",
    "#     coords: Pandas dataframe containing all coordinates from all frames.\n",
    "#     frame: Frame number for which we're getting coordinates.\n",
    "#\n",
    "# Returns:\n",
    "#     A list of indices corresponding to the coordinates for the specified frame\n",
    "#\n",
    "\n",
    "def get_current_coords(coords, frame):\n",
    "    ### CODE GOES HERE ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### TEST SPACE ### \n",
    "\n",
    "# Loading coordinates\n",
    "path = \"./assets/UntrackedCoordinates.csv\"\n",
    "coords = util.load_coordinates(path);\n",
    "\n",
    "# Testing on coordinates from frame 0\n",
    "frame = 0\n",
    "rows_0 = get_current_coords(coords,frame)\n",
    "print(\"Row indices for points in frame 0:\\n%a\\n\" % rows_0)\n",
    "\n",
    "# Testing on coordinates from frame 1\n",
    "frame = 1\n",
    "rows_1 = get_current_coords(coords,frame)\n",
    "print(\"Row indices for points in frame 1:\\n%a\\n\" % rows_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3 - Assign tracks IDs to coordinates in first frame\n",
    "Eventually, all points will be assigned track IDs.  As we work through all the frames, points will be linked to those in a previous frame and inherit their track IDs.  This way, track ID gets propagated through the timeseries.  To kick things off, we need to assign all points in the first frame unique IDs.  The function we create here will identify any points in a specific frame (for now, frame 0) that don't have assigned track IDs and assign them the smallest, currently unused track ID.  By keeping the function general like this we can reuse it later on when find any points in any frame that didn't get linked back to an existing track.\n",
    "\n",
    "### Aims\n",
    "- Assign unique track ID numbers to the rows in the coordinates dataframe (\"coords\") that match the provided rows\n",
    "- Arguments: The coordinates dataframe and the row indices to assign track IDs for\n",
    "- Return: Nothing!  This function will update the input dataframe\n",
    "\n",
    "### Notes\n",
    "- This function will update the coordinates dataframe that is passed to the function, so there are no returned values.\n",
    "- With Pandas you can get the maximum value in a column using the function \"max()\".  For example, the following will give the maximum x-value in the whole dataframe:\n",
    "\n",
    "```Python\n",
    "    max_x_value = coords.X.max()\n",
    "```\n",
    "\n",
    "- One way to do this is to loop over all specified rows, checking if the current track ID is \"0\".  If it is, give it the next available track ID."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to assign unique track ID numbers to any unassigned coordinates at the specified row\n",
    "# indices.  Unassigned coordinates are identified by a track ID of 0.\n",
    "#\n",
    "# Args:\n",
    "#     coords: Pandas dataframe containing all coordinates from all frames.\n",
    "#     point_rows: list of row indices for tracks to have track IDs assigned.  All, some or \n",
    "#                 none of these may need track IDs assigning.\n",
    "#\n",
    "# Returns:\n",
    "#     This function updates the provided coordinate dataframe, so does not return anything.\n",
    "#\n",
    "\n",
    "def assign_new_IDs(coords, point_rows):\n",
    "    ### CODE GOES HERE ###        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### TEST SPACE ###\n",
    "\n",
    "# Loading coordinates\n",
    "path = \"./assets/UntrackedCoordinates.csv\"\n",
    "coords = util.load_coordinates(path);\n",
    "\n",
    "# Getting all points in the first frame\n",
    "frame = 0\n",
    "rows_0 = get_current_coords(coords,frame)\n",
    "\n",
    "# Assigning unique track IDs to all these points\n",
    "assign_new_IDs(coords,rows_0)\n",
    "\n",
    "# Displaying the first 30 rows of our coordinate dataframe\n",
    "print(coords.head(30))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4a - Getting all tracks\n",
    "When assigning track links for each frame, we need to know what tracks are available to link to.  Here, we want to identify the rows corresponding to the most recent instance of each track.\n",
    "\n",
    "To start with, we're going to consider any track that has existed as being \"available\".  As such, we just want to get the most recent row of the coordinates dataframe for each track.  In the second half of this exercise (4b) we will introduce a limit to how far back we will search for the most recent instance of each track.\n",
    "\n",
    "### Aims\n",
    "- Get a list containing the row indices of the most recent coordinate for each track\n",
    "- Arguments: The coordinates dataframe\n",
    "- Return: A list of row indices for the most recent instance of each track\n",
    "\n",
    "### Notes\n",
    "- To check our code works it's useful to have some coordinates that have already been tracked.  Therefore, for this exercise, the \"coords\" dataframe we load has been tracked up to frame 49\n",
    "- We will need to find the unique track IDs in our coordinates dataframe.  It's possible to get all the unique values in a list using Numpy's \"unique\" function.  For example, to get the unique values in the list \"my_list\" we would call:\n",
    "```Python\n",
    "    unique_values = np.unique(my_list)\n",
    "```\n",
    "- You can get the last element of a list using *-1:* as the index.  For example, the following will create a list containing three numbers, then print the value of the final row (i.e. \"21\"):\n",
    "\n",
    "```Python\n",
    "    demo_array = np.array([42,63,21])\n",
    "    print(\"Final element = \",a[-1:])\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to find the row indices corresponding to the most recent instance of each track.\n",
    "#\n",
    "# Args:\n",
    "#     coords: Pandas dataframe containing all coordinates from all frames.\n",
    "#\n",
    "# Returns:\n",
    "#     A list containing the row indices for the most recent instance of each track.\n",
    "#\n",
    "\n",
    "def get_all_tracks(coords):    \n",
    "    ### CODE GOES HERE ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### TEST SPACE ### \n",
    "\n",
    "# Loading coordinates\n",
    "path = \"./assets/PartiallyTrackedCoordinates.csv\"\n",
    "coords = util.load_coordinates(path);\n",
    "\n",
    "# Getting row indices for all available tracks\n",
    "track_rows = get_all_tracks(coords)\n",
    "\n",
    "# Displaying the points we have identified\n",
    "print(\"Identified %i available tracks\" % len(track_rows))\n",
    "print(\"\")\n",
    "print(coords.loc[track_rows])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4b - Getting all RECENT  tracks\n",
    "In the first half of this exercise we identified the most recent instance of each track.  However, we may also choose to only allow links to tracks identified within a specific number of frames.  For example, we may not want to allow a point in frame 42 to link back to a track last seen in frame 9.\n",
    "\n",
    "### Aims\n",
    "- Repeat the code from Exercise 4a, but only permit links to tracks between the provided values of *start_frame* and *end_frame*\n",
    "- Arguments: The coordinates dataframe, the start and end frames, between which we will get available tracks\n",
    "- Return: A list of row indices for the most recent instance of each track\n",
    "\n",
    "### Notes\n",
    "- As before, we'll load the partially-tracked coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to find the row indices corresponding to the most recent instance of each track.\n",
    "# This will only return row indices for tracks present within a specified frame range.\n",
    "#\n",
    "# Args:\n",
    "#     coords: Pandas dataframe containing all coordinates from all frames\n",
    "#     start_frame: First frame of range for which tracks are considered \"available\" for linking to\n",
    "#     end_frame: Final frame of range for which tracks are considered \"available\" for linking to.  This\n",
    "#                will typically be the frame immediately prior to the \"current\" frame.\n",
    "#\n",
    "# Returns:\n",
    "#     A list containing the row indices for the most recent instance of each available track.\n",
    "#\n",
    "\n",
    "def get_available_tracks(coords, start_frame, end_frame):\n",
    "    ### CODE GOES HERE ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### TEST SPACE ### \n",
    "\n",
    "# Loading coordinates\n",
    "path = \"./assets/PartiallyTrackedCoordinates.csv\"\n",
    "coords = util.load_coordinates(path);\n",
    "\n",
    "# Getting row indices for all available tracks\n",
    "start_frame = 44\n",
    "end_frame = 49\n",
    "track_rows = get_available_tracks(coords,start_frame,end_frame)\n",
    "\n",
    "# Displaying the points we have identified\n",
    "print(\"Identified %i available tracks\" % len(track_rows))\n",
    "print(\"\")\n",
    "print(coords.loc[track_rows])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 5 - Calculating cost for linking two points\n",
    "In a couple of steps we will end up with two lists: one for points in the current frame; the other for the most recent point in all available tracks.  Links between the points in the two lists will be assigned based on the cost associated with making that link.  Here, we want to create a function that will calculate the cost associated with two given points.  To start with, we will just calculate the cost as the distance between the two points, but we could also add costs associated with other metrics.  For example, we may want to penalise links which would see the size or intensity of the object change too much.\n",
    "\n",
    "### Aims\n",
    "- Calculate the cost for linking two provided coordinates\n",
    "- Arguments: The coordinates dataframe, the row index for last instance of the track, row index for current point, a distance threshold\n",
    "- Return: A single, floating point value corresponding to the cost for linking these two points\n",
    "\n",
    "### Notes\n",
    "- For this we'll load in the full set of coordinates, but only calculate the cost for the first points in frame 1 and frame 2.\n",
    "- In this case, the \"cost\" is simply the straight-line distance between the two points.\n",
    "- To get the distance between two points we can use the following equation, where x1,y1 and x2,y2 are the xy coordinates of the two points:\n",
    "\n",
    "```Python\n",
    "    dx = x1-x2\n",
    "    dy = y1-y2\n",
    "    dist = math.sqrt(dx*dx + dy*dy)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to calculate the cost of linking two points.  This function simply calculates the cost as \n",
    "# the distance between the two points.\n",
    "#\n",
    "# Args:\n",
    "#     coords: Pandas dataframe containing all coordinates from all frames\n",
    "#     track_row: Row index for the most recent coordinate in a track\n",
    "#     point_row: Row index for a coordinate in the current frame\n",
    "#     thresh: Maximum permitted distance between two points.  If points are separated by more than this\n",
    "#             the cost will be assigned a very large value to ensure a very low probability of the link\n",
    "#             being assigned by the Munkres algorithm.\n",
    "#\n",
    "# Returns:\n",
    "#     A floating point (decimal) value corresponding to the cost of linking the two specified points.\n",
    "#\n",
    "\n",
    "def calculate_cost(coords,track_row,point_row,thresh):\n",
    "    ### CODE GOES HERE ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### TEST SPACE ### \n",
    "\n",
    "# Loading coordinates\n",
    "path = \"./assets/UntrackedCoordinates.csv\"\n",
    "coords = util.load_coordinates(path);\n",
    "\n",
    "# Let's get the indices for the coordinates in the first frame\n",
    "rows_0 = get_current_coords(coords,0)\n",
    "\n",
    "# For this cost calculation we'll just use the first coordinate in this frame\n",
    "row_0 = rows_0[0]\n",
    "\n",
    "# Let's take a look at that coordinate\n",
    "print(\"Coordinate 1:\\n%a\\n\" % coords.loc[row_0])\n",
    "\n",
    "# Now, let's do the same for the second frame\n",
    "rows_1 = get_current_coords(coords,1)\n",
    "row_1 = rows_1[0]\n",
    "print(\"Coordinate 2:\\n%a\\n\" % coords.loc[row_1])\n",
    "\n",
    "# Calculating the cost of the two points with the threshold set high\n",
    "thresh = 10;\n",
    "cost = calculate_cost(coords,row_0,row_1,thresh)\n",
    "print(\"Cost = %f (thresh = %f)\" % (cost,thresh))\n",
    "\n",
    "# Now, we'll do the same, but with a lower threshold\n",
    "thresh = 0.5;\n",
    "cost = calculate_cost(coords,row_0,row_1,thresh)\n",
    "print(\"Cost = %f (thresh = %f)\" % (cost,thresh))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 6 - Calculate cost matrix\n",
    "Now we can calculate the cost for a single pair of points we need to do it for all point pairs.  We will create a function which takes the points in the current frame and the available track points, then generates a 2D cost matrix.  The cost matrix will have a column for each point in the current frame and a row for each available track.  The value of each element will therefore be the cost of linking the point and track corresponding to that column and row.  We also want to limit the distance that links can be made over; therefore, any point-track pairs separations that are greater than a specific distance will be set to infinity.  This doesn't prevent them being linked (if there are no better options, the assignment algorithm will still suggest that as a link), but it makes it less likely.\n",
    "\n",
    "Once we have the cost matrix, we can use SciPy's linear_sum_assignment to calculate the optimal assignments.  This step has been included as the final step in the test space.  We will look into the format of the assignment output in more detail as part of exercise 7.\n",
    "\n",
    "### Aims\n",
    "- Given two lists containing row indices, the first for the most recent instance of all available tracks and the second for all points in the current frame, calculate a cost matrix.\n",
    "- Arguments: The coordinates dataframe, the row indices for tracks, row indices for current points, a distance threshold\n",
    "- Return: An array of floating point values corresponding to the costs for linking all points\n",
    "\n",
    "### Notes\n",
    "- The cost matrix will need to be a 2D Numpy array with a row for every available track and a column for every point in the current frame.\n",
    "- One way to do this is to use nested loops.  \"Nesting loops\" simply means to put one loop inside another.  The following example has two nested loops, which print their iteration numbers:\n",
    "\n",
    "```Python\n",
    "    for i in range(0,3):\n",
    "        for j in range(7,9):\n",
    "            print(i,\",\",j)\n",
    "``` \n",
    "        This will print the following: \"0,7\", \"0,8\", \"1,7\", \"1,8\", \"2,7\", \"2,8\"\n",
    "    \n",
    "- When looping over a structure you can use the *enumerate* function get the current iteration index as well as the current value:\n",
    "\n",
    "```Python\n",
    "    demo_array = np.array([42,63,21])\n",
    "    for i,val in enumerate(demo_array):\n",
    "        print(\"Iter=%i,value=%i\" % (i,val))\n",
    "```\n",
    "        This will print \"Iter=0,value=42\", \"Iter=1,value=63\", \"Iter=2,value=21\"\n",
    "\n",
    "- So we can easily see how this is working, we'll load a special set of coordinates.  This set only has 3 points in the first frame and 4 points in the second frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to calculate the cost of linking all points specified by two lists of row indices.  The\n",
    "# first list of indices corresponds to the most recent coordinates in each available track.  The \n",
    "# second list of indices corresponds to the coordinates in the current frame.  Costs are added to a \n",
    "# 2D Numpy array.\n",
    "#\n",
    "# Args:\n",
    "#     coords: Pandas dataframe containing all coordinates from all frames\n",
    "#     track_rows: Row indices for the most recent coordinate in all available tracks\n",
    "#     point_rows: Row indices for all coordinates in the current frame\n",
    "#     thresh: Maximum permitted distance between two points.  If points are separated by more than this\n",
    "#             the cost will be assigned a very large value to ensure a very low probability of the link\n",
    "#             being assigned by the Munkres algorithm.\n",
    "#\n",
    "# Returns:\n",
    "#     A 2D Numpy array containing all costs.  Each row of the array corresponds to a track and each column\n",
    "#     to a current coordinate.  The intersection of each row and column is the associated cost.\n",
    "#\n",
    "\n",
    "def calculate_cost_matrix(coords,track_rows,point_rows,thresh):\n",
    "    ### CODE GOES HERE ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### TEST SPACE ### \n",
    "\n",
    "# Loading coordinates\n",
    "path = \"./assets/TestCoordinatesForCostMatrix.csv\"\n",
    "coords = util.load_coordinates(path);\n",
    "\n",
    "# Let's get the indices for the coordinates in the first (available tracks) and second (current points) frames\n",
    "track_rows = get_current_coords(coords,0)\n",
    "point_rows = get_current_coords(coords,1)\n",
    "\n",
    "# Calculating the cost matrix\n",
    "thresh = 20\n",
    "costs = calculate_cost_matrix(coords,track_rows,point_rows,thresh)\n",
    "\n",
    "# Displaying our costs\n",
    "print(\"Cost matrix:\\n%a\\n\" % costs)\n",
    "\n",
    "# Using SciPy's linear_sum_assignment solver to get the optimal assignments\n",
    "assignments = linear_sum_assignment(costs)\n",
    "print(\"Assignments:\\n%a\\n%a\\n\" % (assignments[0],assignments[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 7 - Assigning tracks\n",
    "As we saw at the end of the last exercise (in the test space), the job of calculating the assignments is done using the Munkres (aka. Kuhn-Munkres or Hungarian) algorithm.  Rather than create our own implementation of this, we can use SciPy's linear_sum_assignment function.  This algorithm simply takes the cost matrix we just created and outputs a Nx2 array of assignments, where N is the number of assignments.  For each row of the assignment result the first column value corresponds to the track and the second column to the linked coordinate.  \n",
    "\n",
    "As mentioned previously, the Munkres algorithm will still assign links that we set to \"infinity\", so before we copy over any track IDs we want to double check the cost is not infinity.\n",
    "\n",
    "### Aims\n",
    "- Using the calculated track-point assignments, apply the relevant track IDs to each point.\n",
    "- Arguments: The coordinates dataframe, the row indices for the tracks and points to be (potentially) linked, the cost matrix for these tracks and points.\n",
    "- Return: Nothing!  This function will update the input dataframe\n",
    "\n",
    "### Notes\n",
    "- If we were to iterate over the *assignments* array we would receive the first column, then second.  Instead, we want to get the rows one at a time.  We can do this using the *zip* function.\n",
    "- For this example we'll go back to using the partially-tracked example coordinates.  This dataset has tracks assigned up to frame 49.\n",
    "- The test space for this exercise uses many of the functions we've created so far: \n",
    "    - Getting available tracks\n",
    "    - Getting coordinates from the current frame\n",
    "    - Calculating the cost matrix\n",
    "    - Calculating assignments\n",
    "    - Inheriting track IDs from assigned tracks\n",
    "    - Creating new track IDs for unlinked points\n",
    "- We'll only permit linking to tracks last seen within 3 frames of the current frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to assign links between two sets of coordinates based on provided costs.\n",
    "#\n",
    "# Args:\n",
    "#     coords: Pandas dataframe containing all coordinates from all frames\n",
    "#     track_rows: Row indices for the most recent coordinate in all available tracks\n",
    "#     point_rows: Row indices for all coordinates in the current frame\n",
    "#     costs: 2D Numpy array containing all costs.  Each row of the array corresponds to a track and each \n",
    "#            column to a current coordinate.  The intersection of each row and column is the associated cost.\n",
    "#\n",
    "# Returns:\n",
    "#     This function updates the provided coordinate dataframe, so does not return anything.\n",
    "#\n",
    "\n",
    "def assign_IDs(coords, track_rows, point_rows, costs):\n",
    "    ### CODE GOES HERE ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### TEST SPACE ### \n",
    "\n",
    "# Loading coordinates\n",
    "path = \"./assets/PartiallyTrackedCoordinates.csv\"\n",
    "coords = util.load_coordinates(path);\n",
    "\n",
    "# Setting some parameters\n",
    "frame = 50\n",
    "frame_thresh = 2\n",
    "linking_thresh = 5\n",
    "\n",
    "# Getting available tracks at frame 11.  We'll only allow links back 5 frames\n",
    "start_frame = frame - frame_thresh\n",
    "end_frame = frame - 1\n",
    "track_rows = get_available_tracks(coords,start_frame,end_frame)\n",
    "\n",
    "# Getting row indices for points in frame 11\n",
    "point_rows = get_current_coords(coords,frame)\n",
    "\n",
    "# Calculating cost matrix\n",
    "costs = calculate_cost_matrix(coords,track_rows,point_rows,linking_thresh)\n",
    "\n",
    "# Calculating assignments using SciPy's linear_sum_assignment\n",
    "assign_IDs(coords, track_rows, point_rows, costs)\n",
    "\n",
    "# Displaying the coordinates for the current frame (with new assigned track IDs)\n",
    "print(\"Coordinates after linking:\\n%a\\n\" % coords.loc[point_rows])\n",
    "\n",
    "# Assigning new track IDs to points that weren't linked\n",
    "assign_new_IDs(coords,point_rows)\n",
    "print(\"Coordinates after new track creation:\\n%a\\n\" % coords.loc[point_rows])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 8 - Putting it all together\n",
    "We now have all the components necessary to construct a full tracking workflow.  For this, we put most of the components in a single for-loop, which will iterate over each frame.  At the end of this we will re-render the overlay to see if our tracking has worked correctly.  Unlike previous exercises, you don't need to create new functions here - you have everything you need.  The aim is to load the full coordinate set, then iterate over each frame, tracking the coordinates.\n",
    "\n",
    "### Aims\n",
    "- Take the various functions you've created over the past few exercises and put them all together, so we can track over all frames.\n",
    "- At the end of this section you want to end up with the *TRACK_ID* column of the *coords* array completed (i.e. all rows have a relevant ID).\n",
    "\n",
    "### Notes\n",
    "- Unlike previous exercises, there are no extra functions to create here.  As such, you'll be adding code direct into the test space.\n",
    "- You may want to look at the test space for exercise 7 as an example of tracking between two frames.\n",
    "- There's a final box of code after the test space which will run the pre-prepared plotting tool.  If the tracking has worked, you should see tracked lines over the images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### TEST SPACE ### \n",
    "\n",
    "# First, a bit of housekeeping, so we can display the overlaid tracks later on\n",
    "%matplotlib widget\n",
    "\n",
    "# Setting parameters\n",
    "np.set_printoptions(precision=3,threshold=sys.maxsize)\n",
    "\n",
    "linking_thresh = 5\n",
    "frame_thresh = 3\n",
    "\n",
    "# Loading image stack\n",
    "path = \"./assets/images/Fluo-N2DH-GOWT1/01/\"\n",
    "images = util.load_images(path);\n",
    "\n",
    "# Loading coordinates\n",
    "path = \"./assets/UntrackedCoordinates.csv\"\n",
    "coords = util.load_coordinates(path);\n",
    "\n",
    "### CODE GOES HERE ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Adding track renders\n",
    "util.show_overlay(images,coords,True)"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}