[Methods]understanding_DLC_pupil_tracking.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "5c0bb549",
   "metadata": {},
   "outputs": [],
   "source": [
    "import datajoint as dj\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from pipeline.pupil import ConfigDeeplabcut,Tracking\n",
    "import os\n",
    "import cv2\n",
    "import deeplabcut as dlc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ffe08168",
   "metadata": {},
   "outputs": [],
   "source": [
    "from deeplabcut.utils import auxiliaryfunctions\n",
    "from pipeline.utils import DLC_tools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "id": "e6daa38a",
   "metadata": {},
   "outputs": [],
   "source": [
    "key = {'animal_id': 26645, 'session': 2, 'scan_idx': 18}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "1e288841",
   "metadata": {},
   "outputs": [],
   "source": [
    "# change config_path if we were to update DLC model configuration\n",
    "temp_config = (ConfigDeeplabcut & dict(\n",
    "    config_path='/mnt/lab/DeepLabCut/pupil_track-Donnie-2019-02-12/config.yaml')).fetch1()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "a9751cca",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save config_path into the key\n",
    "key['config_path'] = temp_config['config_path']\n",
    "\n",
    "config = auxiliaryfunctions.read_config(temp_config['config_path'])\n",
    "config['config_path'] = temp_config['config_path']\n",
    "config['shuffle'] = temp_config['shuffle']\n",
    "config['trainingsetindex'] = temp_config['trainingsetindex']\n",
    "\n",
    "trainFraction = config['TrainingFraction'][config['trainingsetindex']]\n",
    "DLCscorer = auxiliaryfunctions.GetScorerName(\n",
    "    config, config['shuffle'], trainFraction)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "37c96d95",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make needed directories\n",
    "# tracking_dir, _ = self.create_tracking_directory(key)\n",
    "vid_path = Tracking.Deeplabcut().get_video_path(key)\n",
    "vid_dir = os.path.dirname(os.path.normpath(vid_path))\n",
    "            \n",
    "tracking_dir_name = os.path.basename(\n",
    "                os.path.normpath(vid_path)).split('.')[0] + '_tracking'\n",
    "tracking_dir = os.path.join(vid_dir, tracking_dir_name)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "d336e14c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a short video (5 seconds long)\n",
    "# short_video_path, original_width, original_height, mid_frame_num = DLC_tools.make_short_video(\n",
    "#     tracking_dir)\n",
    "\n",
    "case = os.path.basename(os.path.normpath(tracking_dir)).split('_tracking')[0]\n",
    "input_video_path = os.path.join(tracking_dir, case + '.avi')\n",
    "cap = cv2.VideoCapture(input_video_path)\n",
    "\n",
    "original_width = cap.get(cv2.CAP_PROP_FRAME_WIDTH)\n",
    "original_height = cap.get(cv2.CAP_PROP_FRAME_HEIGHT)\n",
    "mid_frame_num = int(cap.get(cv2.CAP_PROP_FRAME_COUNT)/2)\n",
    "\n",
    "suffix = '_short.avi'\n",
    "short_video_path = os.path.join(tracking_dir, 'short', case + suffix)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "d575c6d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# add original width and height to config\n",
    "config['original_width'] = original_width\n",
    "config['original_height'] = original_height\n",
    "\n",
    "# save info about short video\n",
    "key['short_vid_starting_index'] = mid_frame_num\n",
    "\n",
    "short_h5_path = short_video_path.split('.')[0] + DLCscorer + '.h5'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "66b3188d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/mnt/scratch10/Two-Photon/Kayla/2021-12-14_14-28-04/26645_2_00018_beh_tracking/short/26645_2_00018_beh_shortDeepCut_resnet50_pupil_trackFeb12shuffle1_600000.h5'"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "short_h5_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94bf9ef5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# predict using the short video\n",
    "# DLC_tools.predict_labels(short_video_path, config)\n",
    "# gputouse=0\n",
    "# destfolder = os.path.dirname(short_video_path)\n",
    "# dlc.analyze_videos(config=config['config_path'], videos=[short_video_path], videotype='avi', shuffle=config['shuffle'],\n",
    "#                        trainingsetindex=config['trainingsetindex'], gputouse=gputouse, save_as_csv=False, destfolder=destfolder)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "5cd3e987",
   "metadata": {},
   "outputs": [],
   "source": [
    "# obtain the cropping coordinates from the prediciton on short video\n",
    "# this won't work after the fact because it removes this directory after insert\n",
    "# cropped_coords = DLC_tools.obtain_cropping_coords(short_h5_path, DLCscorer, config)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "690764e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pipeline import pupil"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "id": "86d3b5e7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'animal_id': 26645,\n",
       " 'session': 2,\n",
       " 'scan_idx': 18,\n",
       " 'config_path': '/mnt/lab/DeepLabCut/pupil_track-Donnie-2019-02-12/config.yaml',\n",
       " 'short_vid_starting_index': 0}"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "key"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "id": "d7bd7250",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        \n",
       "        <style type=\"text/css\">\n",
       "            .Relation{\n",
       "                border-collapse:collapse;\n",
       "            }\n",
       "            .Relation th{\n",
       "                background: #A0A0A0; color: #ffffff; padding:4px; border:#f0e0e0 1px solid;\n",
       "                font-weight: normal; font-family: monospace; font-size: 100%;\n",
       "            }\n",
       "            .Relation td{\n",
       "                padding:4px; border:#f0e0e0 1px solid; font-size:100%;\n",
       "            }\n",
       "            .Relation tr:nth-child(odd){\n",
       "                background: #ffffff;\n",
       "            }\n",
       "            .Relation tr:nth-child(even){\n",
       "                background: #f3f1ff;\n",
       "            }\n",
       "            /* Tooltip container */\n",
       "            .djtooltip {\n",
       "            }\n",
       "            /* Tooltip text */\n",
       "            .djtooltip .djtooltiptext {\n",
       "                visibility: hidden;\n",
       "                width: 120px;\n",
       "                background-color: black;\n",
       "                color: #fff;\n",
       "                text-align: center;\n",
       "                padding: 5px 0;\n",
       "                border-radius: 6px;\n",
       "                /* Position the tooltip text - see examples below! */\n",
       "                position: absolute;\n",
       "                z-index: 1;\n",
       "            }\n",
       "            #primary {\n",
       "                font-weight: bold;\n",
       "                color: black;\n",
       "            }\n",
       "\n",
       "            #nonprimary {\n",
       "                font-weight: normal;\n",
       "                color: white;\n",
       "            }\n",
       "\n",
       "            /* Show the tooltip text when you mouse over the tooltip container */\n",
       "            .djtooltip:hover .djtooltiptext {\n",
       "                visibility: visible;\n",
       "            }\n",
       "        </style>\n",
       "        \n",
       "        <b></b>\n",
       "            <div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
       "            <table border=\"1\" class=\"Relation\">\n",
       "                <thead> <tr style=\"text-align: right;\"> <th> <div class=\"djtooltip\">\n",
       "                                <p id=\"primary\">animal_id</p>\n",
       "                                <span class=\"djtooltiptext\">id number</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"primary\">session</p>\n",
       "                                <span class=\"djtooltiptext\">session index for the mouse</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"primary\">scan_idx</p>\n",
       "                                <span class=\"djtooltiptext\">number of TIFF stack file</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"primary\">tracking_method</p>\n",
       "                                <span class=\"djtooltiptext\">method used for pupil tracking</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"primary\">frame_id</p>\n",
       "                                <span class=\"djtooltiptext\">frame id with matlab based 1 indexing</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"nonprimary\">center</p>\n",
       "                                <span class=\"djtooltiptext\">center of the circle in (x, y) of image</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"nonprimary\">radius</p>\n",
       "                                <span class=\"djtooltiptext\">radius of the circle</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"nonprimary\">visible_portion</p>\n",
       "                                <span class=\"djtooltiptext\">portion of visible pupil area given a fitted circle frame. Please refer DLC_tools.PupilFitting.detect_visible_pupil_area for more details</span>\n",
       "                            </div> </th> </tr> </thead>\n",
       "                <tbody> <tr> <td>26645</td>\n",
       "<td>2</td>\n",
       "<td>18</td>\n",
       "<td>2</td>\n",
       "<td>0</td>\n",
       "<td>=BLOB=</td>\n",
       "<td>nan</td>\n",
       "<td>-3.0</td></tr><tr><td>26645</td>\n",
       "<td>2</td>\n",
       "<td>18</td>\n",
       "<td>2</td>\n",
       "<td>1</td>\n",
       "<td>=BLOB=</td>\n",
       "<td>nan</td>\n",
       "<td>-3.0</td></tr><tr><td>26645</td>\n",
       "<td>2</td>\n",
       "<td>18</td>\n",
       "<td>2</td>\n",
       "<td>2</td>\n",
       "<td>=BLOB=</td>\n",
       "<td>nan</td>\n",
       "<td>-3.0</td></tr><tr><td>26645</td>\n",
       "<td>2</td>\n",
       "<td>18</td>\n",
       "<td>2</td>\n",
       "<td>3</td>\n",
       "<td>=BLOB=</td>\n",
       "<td>nan</td>\n",
       "<td>-3.0</td></tr><tr><td>26645</td>\n",
       "<td>2</td>\n",
       "<td>18</td>\n",
       "<td>2</td>\n",
       "<td>4</td>\n",
       "<td>=BLOB=</td>\n",
       "<td>nan</td>\n",
       "<td>-3.0</td></tr><tr><td>26645</td>\n",
       "<td>2</td>\n",
       "<td>18</td>\n",
       "<td>2</td>\n",
       "<td>5</td>\n",
       "<td>=BLOB=</td>\n",
       "<td>nan</td>\n",
       "<td>-3.0</td></tr><tr><td>26645</td>\n",
       "<td>2</td>\n",
       "<td>18</td>\n",
       "<td>2</td>\n",
       "<td>6</td>\n",
       "<td>=BLOB=</td>\n",
       "<td>nan</td>\n",
       "<td>-3.0</td></tr><tr><td>26645</td>\n",
       "<td>2</td>\n",
       "<td>18</td>\n",
       "<td>2</td>\n",
       "<td>7</td>\n",
       "<td>=BLOB=</td>\n",
       "<td>nan</td>\n",
       "<td>-3.0</td></tr><tr><td>26645</td>\n",
       "<td>2</td>\n",
       "<td>18</td>\n",
       "<td>2</td>\n",
       "<td>8</td>\n",
       "<td>=BLOB=</td>\n",
       "<td>nan</td>\n",
       "<td>-3.0</td></tr><tr><td>26645</td>\n",
       "<td>2</td>\n",
       "<td>18</td>\n",
       "<td>2</td>\n",
       "<td>9</td>\n",
       "<td>=BLOB=</td>\n",
       "<td>nan</td>\n",
       "<td>-3.0</td></tr><tr><td>26645</td>\n",
       "<td>2</td>\n",
       "<td>18</td>\n",
       "<td>2</td>\n",
       "<td>10</td>\n",
       "<td>=BLOB=</td>\n",
       "<td>nan</td>\n",
       "<td>-3.0</td></tr><tr><td>26645</td>\n",
       "<td>2</td>\n",
       "<td>18</td>\n",
       "<td>2</td>\n",
       "<td>11</td>\n",
       "<td>=BLOB=</td>\n",
       "<td>nan</td>\n",
       "<td>-3.0</td> </tr> </tbody>\n",
       "            </table>\n",
       "            <p>...</p>\n",
       "            <p>Total: 148846</p></div>\n",
       "            "
      ],
      "text/plain": [
       "*animal_id    *session    *scan_idx    *tracking_meth *frame_id    center     radius     visible_portio\n",
       "+-----------+ +---------+ +----------+ +------------+ +----------+ +--------+ +--------+ +------------+\n",
       "26645         2           18           2              0            =BLOB=     nan        -3.0          \n",
       "26645         2           18           2              1            =BLOB=     nan        -3.0          \n",
       "26645         2           18           2              2            =BLOB=     nan        -3.0          \n",
       "26645         2           18           2              3            =BLOB=     nan        -3.0          \n",
       "26645         2           18           2              4            =BLOB=     nan        -3.0          \n",
       "26645         2           18           2              5            =BLOB=     nan        -3.0          \n",
       "26645         2           18           2              6            =BLOB=     nan        -3.0          \n",
       "26645         2           18           2              7            =BLOB=     nan        -3.0          \n",
       "26645         2           18           2              8            =BLOB=     nan        -3.0          \n",
       "26645         2           18           2              9            =BLOB=     nan        -3.0          \n",
       "26645         2           18           2              10           =BLOB=     nan        -3.0          \n",
       "26645         2           18           2              11           =BLOB=     nan        -3.0          \n",
       "   ...\n",
       " (Total: 148846)"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pupil.FittedPupil.Circle & key"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "id": "f69fe5aa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        \n",
       "        <style type=\"text/css\">\n",
       "            .Relation{\n",
       "                border-collapse:collapse;\n",
       "            }\n",
       "            .Relation th{\n",
       "                background: #A0A0A0; color: #ffffff; padding:4px; border:#f0e0e0 1px solid;\n",
       "                font-weight: normal; font-family: monospace; font-size: 100%;\n",
       "            }\n",
       "            .Relation td{\n",
       "                padding:4px; border:#f0e0e0 1px solid; font-size:100%;\n",
       "            }\n",
       "            .Relation tr:nth-child(odd){\n",
       "                background: #ffffff;\n",
       "            }\n",
       "            .Relation tr:nth-child(even){\n",
       "                background: #f3f1ff;\n",
       "            }\n",
       "            /* Tooltip container */\n",
       "            .djtooltip {\n",
       "            }\n",
       "            /* Tooltip text */\n",
       "            .djtooltip .djtooltiptext {\n",
       "                visibility: hidden;\n",
       "                width: 120px;\n",
       "                background-color: black;\n",
       "                color: #fff;\n",
       "                text-align: center;\n",
       "                padding: 5px 0;\n",
       "                border-radius: 6px;\n",
       "                /* Position the tooltip text - see examples below! */\n",
       "                position: absolute;\n",
       "                z-index: 1;\n",
       "            }\n",
       "            #primary {\n",
       "                font-weight: bold;\n",
       "                color: black;\n",
       "            }\n",
       "\n",
       "            #nonprimary {\n",
       "                font-weight: normal;\n",
       "                color: white;\n",
       "            }\n",
       "\n",
       "            /* Show the tooltip text when you mouse over the tooltip container */\n",
       "            .djtooltip:hover .djtooltiptext {\n",
       "                visibility: visible;\n",
       "            }\n",
       "        </style>\n",
       "        \n",
       "        <b></b>\n",
       "            <div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
       "            <table border=\"1\" class=\"Relation\">\n",
       "                <thead> <tr style=\"text-align: right;\"> <th> <div class=\"djtooltip\">\n",
       "                                <p id=\"primary\">animal_id</p>\n",
       "                                <span class=\"djtooltiptext\">id number</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"primary\">session</p>\n",
       "                                <span class=\"djtooltiptext\">session index for the mouse</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"primary\">scan_idx</p>\n",
       "                                <span class=\"djtooltiptext\">number of TIFF stack file</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"primary\">tracking_method</p>\n",
       "                                <span class=\"djtooltiptext\">method used for pupil tracking</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"nonprimary\">short_vid_starting_index</p>\n",
       "                                <span class=\"djtooltiptext\">middle frame index of the original video</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"nonprimary\">cropped_x0</p>\n",
       "                                <span class=\"djtooltiptext\">start width coord wrt original video</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"nonprimary\">cropped_x1</p>\n",
       "                                <span class=\"djtooltiptext\">end width coord wrt original video</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"nonprimary\">cropped_y0</p>\n",
       "                                <span class=\"djtooltiptext\">start height coord wrt original video</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"nonprimary\">cropped_y1</p>\n",
       "                                <span class=\"djtooltiptext\">end height coord wrt original video</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"nonprimary\">added_pixels</p>\n",
       "                                <span class=\"djtooltiptext\">number of pixels added around the cropping coords</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"nonprimary\">config_path</p>\n",
       "                                <span class=\"djtooltiptext\">path to deeplabcut config yaml</span>\n",
       "                            </div> </th> </tr> </thead>\n",
       "                <tbody> <tr> <td>26645</td>\n",
       "<td>2</td>\n",
       "<td>18</td>\n",
       "<td>2</td>\n",
       "<td>74423</td>\n",
       "<td>23</td>\n",
       "<td>365</td>\n",
       "<td>0</td>\n",
       "<td>281</td>\n",
       "<td>100</td>\n",
       "<td>/mnt/lab/DeepLabCut/pupil_track-Donnie-2019-02-12/config.yaml</td> </tr> </tbody>\n",
       "            </table>\n",
       "            \n",
       "            <p>Total: 1</p></div>\n",
       "            "
      ],
      "text/plain": [
       "*animal_id    *session    *scan_idx    *tracking_meth short_vid_star cropped_x0     cropped_x1     cropped_y0     cropped_y1     added_pixels   config_path   \n",
       "+-----------+ +---------+ +----------+ +------------+ +------------+ +------------+ +------------+ +------------+ +------------+ +------------+ +------------+\n",
       "26645         2           18           2              74423          23             365            0              281            100            /mnt/lab/DeepL\n",
       " (Total: 1)"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Tracking.Deeplabcut & key"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "id": "ef7fa091",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        \n",
       "        <style type=\"text/css\">\n",
       "            .Relation{\n",
       "                border-collapse:collapse;\n",
       "            }\n",
       "            .Relation th{\n",
       "                background: #A0A0A0; color: #ffffff; padding:4px; border:#f0e0e0 1px solid;\n",
       "                font-weight: normal; font-family: monospace; font-size: 100%;\n",
       "            }\n",
       "            .Relation td{\n",
       "                padding:4px; border:#f0e0e0 1px solid; font-size:100%;\n",
       "            }\n",
       "            .Relation tr:nth-child(odd){\n",
       "                background: #ffffff;\n",
       "            }\n",
       "            .Relation tr:nth-child(even){\n",
       "                background: #f3f1ff;\n",
       "            }\n",
       "            /* Tooltip container */\n",
       "            .djtooltip {\n",
       "            }\n",
       "            /* Tooltip text */\n",
       "            .djtooltip .djtooltiptext {\n",
       "                visibility: hidden;\n",
       "                width: 120px;\n",
       "                background-color: black;\n",
       "                color: #fff;\n",
       "                text-align: center;\n",
       "                padding: 5px 0;\n",
       "                border-radius: 6px;\n",
       "                /* Position the tooltip text - see examples below! */\n",
       "                position: absolute;\n",
       "                z-index: 1;\n",
       "            }\n",
       "            #primary {\n",
       "                font-weight: bold;\n",
       "                color: black;\n",
       "            }\n",
       "\n",
       "            #nonprimary {\n",
       "                font-weight: normal;\n",
       "                color: white;\n",
       "            }\n",
       "\n",
       "            /* Show the tooltip text when you mouse over the tooltip container */\n",
       "            .djtooltip:hover .djtooltiptext {\n",
       "                visibility: visible;\n",
       "            }\n",
       "        </style>\n",
       "        \n",
       "        <b>Minimal info needed to load deeplabcut model</b>\n",
       "            <div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
       "            <table border=\"1\" class=\"Relation\">\n",
       "                <thead> <tr style=\"text-align: right;\"> <th> <div class=\"djtooltip\">\n",
       "                                <p id=\"primary\">config_path</p>\n",
       "                                <span class=\"djtooltiptext\">path to deeplabcut config</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"nonprimary\">shuffle</p>\n",
       "                                <span class=\"djtooltiptext\">shuffle number used for the trained dlc model. Needed for dlc.analyze_videos</span>\n",
       "                            </div></th><th><div class=\"djtooltip\">\n",
       "                                <p id=\"nonprimary\">trainingsetindex</p>\n",
       "                                <span class=\"djtooltiptext\">trainingset index used for the trained dlc. model. Needed for dlc.analyze_videos</span>\n",
       "                            </div> </th> </tr> </thead>\n",
       "                <tbody> <tr> <td>/mnt/lab/DeepLabCut/pupil_track-Donnie-2019-02-12/config.yaml</td>\n",
       "<td>1</td>\n",
       "<td>0</td> </tr> </tbody>\n",
       "            </table>\n",
       "            \n",
       "            <p>Total: 1</p></div>\n",
       "            "
      ],
      "text/plain": [
       "*config_path   shuffle     trainingsetind\n",
       "+------------+ +---------+ +------------+\n",
       "/mnt/lab/DeepL 1           0             \n",
       " (Total: 1)"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ConfigDeeplabcut & (Tracking.Deeplabcut & key)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "0b39711d",
   "metadata": {},
   "outputs": [],
   "source": [
    "dlc_config = (ConfigDeeplabcut & (\n",
    "                Tracking.Deeplabcut & key)).fetch1()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "d09aa8a1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'config_path': '/mnt/lab/DeepLabCut/pupil_track-Donnie-2019-02-12/config.yaml',\n",
       " 'shuffle': 1,\n",
       " 'trainingsetindex': 0}"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dlc_config"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "id": "6006208a",
   "metadata": {},
   "outputs": [],
   "source": [
    "config = auxiliaryfunctions.read_config(dlc_config['config_path'])\n",
    "config['config_path'] = dlc_config['config_path']\n",
    "config['shuffle'] = dlc_config['shuffle']\n",
    "config['trainingsetindex'] = dlc_config['trainingsetindex']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "id": "bc222e96",
   "metadata": {},
   "outputs": [],
   "source": [
    "avi_path = (pupil.Eye & key).get_video_path()\n",
    "nframes = (pupil.Eye & key).fetch1('total_frames')\n",
    "\n",
    "# find path to original video symlink\n",
    "base_path = os.path.splitext(avi_path)[0] + '_tracking'\n",
    "video_path = os.path.join(base_path, os.path.basename(avi_path))\n",
    "\n",
    "config['orig_video_path'] = video_path\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "id": "a169ad1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# find croppoing coords\n",
    "cropped_coords = (Tracking.Deeplabcut & key).fetch1(\n",
    "    'cropped_x0', 'cropped_x1', 'cropped_y0', 'cropped_y1')\n",
    "\n",
    "config['cropped_coords'] = cropped_coords\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69304eed",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b7b7061",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "pupil_fit = DLC_tools.DeeplabcutPupilFitting(\n",
    "    config=config, bodyparts='all', cropped=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b40c2f48",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15aed962",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4d536ef",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bba8d430",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6215e954",
   "metadata": {},
   "outputs": [],
   "source": [
    "@schema\n",
    "class Tracking(dj.Computed):\n",
    "    definition = \"\"\"\n",
    "    -> Eye\n",
    "    -> shared.TrackingMethod\n",
    "    ---\n",
    "    tracking_ts=CURRENT_TIMESTAMP   : timestamp  # automatic\n",
    "    \"\"\"\n",
    "\n",
    "    class Deeplabcut(dj.Part):\n",
    "        definition = \"\"\"\n",
    "        -> master\n",
    "        ---\n",
    "        short_vid_starting_index    : int unsigned          # middle frame index of the original video\n",
    "        cropped_x0                  : smallint unsigned     # start width coord wrt original video\n",
    "        cropped_x1                  : smallint unsigned     # end width coord wrt original video\n",
    "        cropped_y0                  : smallint unsigned     # start height coord wrt original video\n",
    "        cropped_y1                  : smallint unsigned     # end height coord wrt original video\n",
    "        added_pixels                : smallint unsigned     # number of pixels added around the cropping coords\n",
    "        config_path                 : varchar(128)          # path to deeplabcut config yaml        \n",
    "        \"\"\"\n",
    "\n",
    "        def get_video_path(self, key):\n",
    "            \"\"\"\n",
    "            Input:\n",
    "                key: dictionary\n",
    "                    A key that consists of animal_id, session, and scan_idx\n",
    "            \"\"\"\n",
    "            video_info = (experiment.Session() *\n",
    "                          experiment.Scan.EyeVideo() & key).fetch1()\n",
    "            video_path = lab.Paths().get_local_path(\n",
    "                \"{behavior_path}/{filename}\".format(**video_info))\n",
    "            return video_path\n",
    "\n",
    "        def create_tracking_directory(self, key):\n",
    "            \"\"\"\n",
    "            this function creates the following directory structure:\n",
    "            video_original_dir\n",
    "                |\n",
    "                |------ video_original\n",
    "                |------ tracking_dir (create_tracking_directory)\n",
    "                            |------- symlink to video_original (add_symlink)\n",
    "                            |------- compressed_cropped_dir\n",
    "                                        |------- cropped_video (generated by make_compressed_cropped_video function)\n",
    "                                        |------- h5 file for cropped video (generated by deeplabcut)\n",
    "                                        |------- pickle for cropped video (generated by deeplabcut)\n",
    "                            |------- short_dir\n",
    "                                        |------- short_video (generated by make_short_video function)\n",
    "                                        |------- h5 file for short video(generated by deeplabcut)\n",
    "                                        |------- pickle for short video (generated by deeplabcut)\n",
    "            Input:\n",
    "                key: dictionary\n",
    "                    a dictionary that contains mouse id, session, and scan idx.\n",
    "            Return:\n",
    "                tracking_dir: string\n",
    "                    a string that specifies the path to the tracking directory\n",
    "            \"\"\"\n",
    "\n",
    "            print(\"Generating tracking directory for \", key)\n",
    "\n",
    "            vid_path = self.get_video_path(key)\n",
    "            vid_dir = os.path.dirname(os.path.normpath(vid_path))\n",
    "            tracking_dir_name = os.path.basename(\n",
    "                os.path.normpath(vid_path)).split('.')[0] + '_tracking'\n",
    "\n",
    "            tracking_dir = os.path.join(vid_dir, tracking_dir_name)\n",
    "\n",
    "            hardlink_path = os.path.join(\n",
    "                tracking_dir, os.path.basename(os.path.normpath(vid_path)))\n",
    "\n",
    "            if not os.path.exists(tracking_dir):\n",
    "\n",
    "                os.mkdir(tracking_dir)\n",
    "                os.mkdir(os.path.join(tracking_dir, 'compressed_cropped'))\n",
    "                os.mkdir(os.path.join(tracking_dir, 'short'))\n",
    "\n",
    "                os.link(vid_path, hardlink_path)\n",
    "\n",
    "            else:\n",
    "                print('{} already exists!'.format(tracking_dir))\n",
    "                print('Removing existing tracking directory and recreating')\n",
    "\n",
    "                shutil.rmtree(tracking_dir)\n",
    "\n",
    "                os.mkdir(tracking_dir)\n",
    "                os.mkdir(os.path.join(tracking_dir, 'compressed_cropped'))\n",
    "                os.mkdir(os.path.join(tracking_dir, 'short'))\n",
    "\n",
    "                os.link(vid_path, hardlink_path)\n",
    "\n",
    "            return tracking_dir, hardlink_path\n",
    "\n",
    "        def make(self, key):\n",
    "            \"\"\"\n",
    "            Use Deeplabcut to label pupil and eyelids\n",
    "            \"\"\"\n",
    "\n",
    "            print('Tracking labels with Deeplabcut!')\n",
    "\n",
    "            # change config_path if we were to update DLC model configuration\n",
    "            temp_config = (ConfigDeeplabcut & dict(\n",
    "                config_path='/mnt/lab/DeepLabCut/pupil_track-Donnie-2019-02-12/config.yaml')).fetch1()\n",
    "\n",
    "            # save config_path into the key\n",
    "            key['config_path'] = temp_config['config_path']\n",
    "\n",
    "            config = auxiliaryfunctions.read_config(temp_config['config_path'])\n",
    "            config['config_path'] = temp_config['config_path']\n",
    "            config['shuffle'] = temp_config['shuffle']\n",
    "            config['trainingsetindex'] = temp_config['trainingsetindex']\n",
    "\n",
    "            trainFraction = config['TrainingFraction'][config['trainingsetindex']]\n",
    "            DLCscorer = auxiliaryfunctions.GetScorerName(\n",
    "                config, config['shuffle'], trainFraction)\n",
    "\n",
    "            # make needed directories\n",
    "            tracking_dir, _ = self.create_tracking_directory(key)\n",
    "\n",
    "            # make a short video (5 seconds long)\n",
    "            short_video_path, original_width, original_height, mid_frame_num = DLC_tools.make_short_video(\n",
    "                tracking_dir)\n",
    "\n",
    "            # add original width and height to config\n",
    "            config['original_width'] = original_width\n",
    "            config['original_height'] = original_height\n",
    "\n",
    "            # save info about short video\n",
    "            key['short_vid_starting_index'] = mid_frame_num\n",
    "\n",
    "            short_h5_path = short_video_path.split('.')[0] + DLCscorer + '.h5'\n",
    "\n",
    "            # predict using the short video\n",
    "            DLC_tools.predict_labels(short_video_path, config)\n",
    "\n",
    "            # obtain the cropping coordinates from the prediciton on short video\n",
    "            cropped_coords = DLC_tools.obtain_cropping_coords(\n",
    "                short_h5_path, DLCscorer, config)\n",
    "\n",
    "            # add 100 pixels around cropping coords. Ensure that it is within the original dim\n",
    "            pixel_num = 100\n",
    "            cropped_coords = DLC_tools.add_pixels(cropped_coords=cropped_coords,\n",
    "                                                  original_width=original_width,\n",
    "                                                  original_height=original_height,\n",
    "                                                  pixel_num=pixel_num)\n",
    "\n",
    "            # make a compressed and cropped video\n",
    "            compressed_cropped_video_path = DLC_tools.make_compressed_cropped_video(\n",
    "                tracking_dir, cropped_coords)\n",
    "\n",
    "            # predict using the compressed and cropped video\n",
    "            DLC_tools.predict_labels(compressed_cropped_video_path, config)\n",
    "\n",
    "            key = dict(key, cropped_x0=cropped_coords['cropped_x0'],\n",
    "                       cropped_x1=cropped_coords['cropped_x1'],\n",
    "                       cropped_y0=cropped_coords['cropped_y0'],\n",
    "                       cropped_y1=cropped_coords['cropped_y1'],\n",
    "                       added_pixels=pixel_num)\n",
    "\n",
    "            self.insert1(key)\n",
    "\n",
    "            # delete short video directory\n",
    "            shutil.rmtree(os.path.dirname(short_video_path))\n",
    "\n",
    "            # delete compressed and cropped video\n",
    "            os.remove(compressed_cropped_video_path)\n",
    "\n",
    "    def make(self, key):\n",
    "        print(\"Tracking for case {}\".format(key))\n",
    "\n",
    "        if key['tracking_method'] == 1:\n",
    "            self.insert1(key)\n",
    "            self.ManualTracking().make(key)\n",
    "        elif key['tracking_method'] == 2:\n",
    "            self.insert1(key)\n",
    "            self.Deeplabcut().make(key)\n",
    "        else:\n",
    "            msg = 'Unrecognized Tracking method {}'.format(\n",
    "                key['tracking_method'])\n",
    "            raise PipelineException(msg)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33a78521",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "    def make(self, key):\n",
    "        print(\"Fitting:\", key)\n",
    "\n",
    "        self.insert1(key)\n",
    "\n",
    "        avi_path = (Eye & key).get_video_path()\n",
    "        nframes = (Eye & key).fetch1('total_frames')\n",
    "\n",
    "        data_circle = []\n",
    "        data_ellipse = []\n",
    "\n",
    "        # deeplabcut 2\n",
    "        elif key['tracking_method'] == 2:\n",
    "\n",
    "            dlc_config = (ConfigDeeplabcut & (\n",
    "                Tracking.Deeplabcut & key)).fetch1()\n",
    "\n",
    "            config = auxiliaryfunctions.read_config(dlc_config['config_path'])\n",
    "            config['config_path'] = dlc_config['config_path']\n",
    "            config['shuffle'] = dlc_config['shuffle']\n",
    "            config['trainingsetindex'] = dlc_config['trainingsetindex']\n",
    "\n",
    "            # find path to original video symlink\n",
    "            base_path = os.path.splitext(avi_path)[0] + '_tracking'\n",
    "            video_path = os.path.join(base_path, os.path.basename(avi_path))\n",
    "\n",
    "            config['orig_video_path'] = video_path\n",
    "\n",
    "            # find croppoing coords\n",
    "            cropped_coords = (Tracking.Deeplabcut & key).fetch1(\n",
    "                'cropped_x0', 'cropped_x1', 'cropped_y0', 'cropped_y1')\n",
    "\n",
    "            config['cropped_coords'] = cropped_coords\n",
    "\n",
    "            pupil_fit = DLC_tools.DeeplabcutPupilFitting(\n",
    "                config=config, bodyparts='all', cropped=True)\n",
    "\n",
    "            for bodypart in pupil_fit.bodyparts:\n",
    "                self.EyePoints.insert1({**key,\n",
    "                    'label': bodypart,\n",
    "                    'x': pupil_fit.df_bodyparts[bodypart]['x'].values,\n",
    "                    'y': pupil_fit.df_bodyparts[bodypart]['y'].values       \n",
    "                })\n",
    "\n",
    "            for frame_num in tqdm(range(nframes)):\n",
    "\n",
    "                fit_dict = pupil_fit.fitted_core(frame_num=frame_num)\n",
    "\n",
    "                # circle info\n",
    "                center = fit_dict['circle_fit']['center']\n",
    "                radius = fit_dict['circle_fit']['radius']\n",
    "                visible_portion = fit_dict['circle_visible']['visible_portion']\n",
    "\n",
    "                data_circle.append([center, radius, visible_portion])\n",
    "\n",
    "                # ellipse info\n",
    "                center = fit_dict['ellipse_fit']['center']\n",
    "                major_radius = fit_dict['ellipse_fit']['major_radius']\n",
    "                minor_radius = fit_dict['ellipse_fit']['minor_radius']\n",
    "                rotation_angle = fit_dict['ellipse_fit']['rotation_angle']\n",
    "                visible_portion = fit_dict['ellipse_visible']['visible_portion']\n",
    "\n",
    "                data_ellipse.append([center, major_radius, minor_radius,\n",
    "                                     rotation_angle, visible_portion])\n",
    "\n",
    "\n",
    "        data_circle = np.array(data_circle)\n",
    "        data_ellipse = np.array(data_ellipse)\n",
    "        # now filter out the outliers by 5.5 std away from mean\n",
    "        rejected_ind = DLC_tools.filter_by_fitting_std(\n",
    "            data=data_circle, fitting_method='circle', std_magnitude=5.5)\n",
    "\n",
    "        data_circle[rejected_ind] = None, None, -3.0\n",
    "\n",
    "        common_entry = np.array(list(key.values()))\n",
    "        common_matrix = np.tile(common_entry, (nframes, 1))\n",
    "\n",
    "        data_circle = np.hstack(\n",
    "            (common_matrix, np.arange(nframes).reshape(-1, 1), data_circle))\n",
    "\n",
    "        # insert data\n",
    "        with dj.config(enable_python_native_blobs = True):\n",
    "            self.Circle.insert(data_circle)\n",
    "\n",
    "        # now repeat the process for ellipse\n",
    "        rejected_ind = DLC_tools.filter_by_fitting_std(\n",
    "            data=data_ellipse, fitting_method='ellipse', std_magnitude=5.5)\n",
    "\n",
    "        data_ellipse[rejected_ind, :] = None, None, None, None, -3.0\n",
    "\n",
    "        data_ellipse = np.hstack(\n",
    "            (common_matrix, np.arange(nframes).reshape(-1, 1), data_ellipse))\n",
    "\n",
    "        with dj.config(enable_python_native_blobs = True):\n",
    "\n",
    "            self.Ellipse.insert(data_ellipse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec950eaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "class DeeplabcutPupilFitting(DeeplabcutPlotBodyparts):\n",
    "    def __init__(self, config, bodyparts='all', cropped=False, filtering=None):\n",
    "        \"\"\"\n",
    "        Input:\n",
    "            config: dictionary\n",
    "                A dictionary that contains animal id, session, scan idx, and a path to config\n",
    "            bodyparts: list\n",
    "                A list that contains bodyparts to plot. Each bodypart is in a string format. If none provided,\n",
    "                then by default it plots ALL existing bodyplots in config.yaml file.\n",
    "        \"\"\"\n",
    "        super().__init__(config, bodyparts=bodyparts, cropped=cropped)\n",
    "\n",
    "        self.complete_eyelid_graph = {'eyelid_top': 'eyelid_top_right',\n",
    "                                      'eyelid_top_right': 'eyelid_right',\n",
    "                                      'eyelid_right': 'eyelid_right_bottom',\n",
    "                                      'eyelid_right_bottom': 'eyelid_bottom',\n",
    "                                      'eyelid_bottom': 'eyelid_bottom_left',\n",
    "                                      'eyelid_bottom_left': 'eyelid_left',\n",
    "                                      'eyelid_left': 'eyelid_left_top',\n",
    "                                      'eyelid_left_top': 'eyelid_top'}\n",
    "\n",
    "        self._circle_threshold_num = 3\n",
    "        self._ellipse_threshold_num = 6\n",
    "        self._circle_color = (0, 255, 0)\n",
    "        self._ellipse_color = (0, 0, 255)\n",
    "\n",
    "    @property\n",
    "    def circle_threshold_num(self):\n",
    "        return self._circle_threshold_num\n",
    "\n",
    "    @circle_threshold_num.setter\n",
    "    def circle_threshold_num(self, value):\n",
    "        if value > len(self.complete_eyelid_graph.keys()):\n",
    "            raise ValueError(\"value must be equal to or less than {}!\".format(\n",
    "                len(self.complete_eyelid_graph.keys())))\n",
    "        else:\n",
    "            self._circle_threshold_num = value\n",
    "\n",
    "    @property\n",
    "    def ellipse_threshold_num(self):\n",
    "        return self._ellipse_threshold_num\n",
    "\n",
    "    @ellipse_threshold_num.setter\n",
    "    def ellipse_threshold_num(self, value):\n",
    "        if value > len(self.complete_eyelid_graph.keys()):\n",
    "            raise ValueError(\"value must be equal to or less than {}!\".format(\n",
    "                len(self.complete_eyelid_graph.keys())))\n",
    "\n",
    "        # 5 is the minimum number needed for ellipse fitting\n",
    "        elif value < 5:\n",
    "            raise ValueError(\"value must be equal to or more than 5!\")\n",
    "        else:\n",
    "            self._ellipse_threshold_num = value\n",
    "\n",
    "    @property\n",
    "    def circle_color(self):\n",
    "        return self._circle_color\n",
    "\n",
    "    @circle_color.setter\n",
    "    def circle_color(self, value):\n",
    "        self._circle_color = value\n",
    "\n",
    "    @property\n",
    "    def ellipse_color(self):\n",
    "        return self._ellipse_color\n",
    "\n",
    "    @ellipse_color.setter\n",
    "    def ellipse_color(self, value):\n",
    "        self._ellipse_color = value\n",
    "\n",
    "    def coords_pcutoff(self, frame_num):\n",
    "        \"\"\"\n",
    "        Given a frame number, return bpindex, x & y coordinates that meet pcutoff criteria\n",
    "        Input:\n",
    "            frame_num: int\n",
    "                A desired frame number\n",
    "        Output:\n",
    "            bpindex: list\n",
    "                A list of integers that match with bodypart. For instance, if the bodypart is ['A','B','C']\n",
    "                and only 'A' and 'C'qualifies the pcutoff, then bpindex = [0,2]\n",
    "            x_coords: pandas series\n",
    "                A pandas series that contains coordinates whose values meet pcutoff criteria\n",
    "            y_coords: pandas series\n",
    "                A pandas series that contains coordinates whose values meet pcutoff criteria\n",
    "        \"\"\"\n",
    "        frame_num_tf = self.tf_likelihood_array[frame_num, :]\n",
    "        bpindex = [i for i, x in enumerate(frame_num_tf) if x]\n",
    "\n",
    "        df_x_coords = self.df_bodyparts_x.loc[frame_num, :][bpindex]\n",
    "        df_y_coords = self.df_bodyparts_y.loc[frame_num, :][bpindex]\n",
    "\n",
    "        return bpindex, df_x_coords, df_y_coords\n",
    "\n",
    "    def connect_eyelids(self, frame_num, frame):\n",
    "        \"\"\"\n",
    "        connect eyelid labels with a straight line. If a label is missing, do not connect and skip to the next label.\n",
    "        Input:\n",
    "            frame_num: int\n",
    "                A desired frame number\n",
    "            frame: numpy array\n",
    "                A frame to be fitted\n",
    "        Output:\n",
    "            A dictionary containing the fitted frame and its corresponding binary mask.\n",
    "            For each key in dictionary:\n",
    "                frame:\n",
    "                    A numpy array frame with eyelids connected\n",
    "                mask:\n",
    "                    A binary numpy array\n",
    "        \"\"\"\n",
    "        mask = np.zeros(frame.shape[:2], dtype=np.uint8)\n",
    "\n",
    "        _, df_x_coords, df_y_coords = self.coords_pcutoff(frame_num)\n",
    "        eyelid_labels = [label for label in list(\n",
    "            df_x_coords.index.get_level_values(0)) if 'eyelid' in label]\n",
    "\n",
    "        for eyelid in eyelid_labels:\n",
    "            next_bp = self.complete_eyelid_graph[eyelid]\n",
    "\n",
    "            if next_bp not in eyelid_labels:\n",
    "                continue\n",
    "\n",
    "            coord_0 = tuple(\n",
    "                map(int, map(round, [df_x_coords[eyelid].values[0], df_y_coords[eyelid].values[0]])))\n",
    "            coord_1 = tuple(\n",
    "                map(int, map(round, [df_x_coords[next_bp].values[0], df_y_coords[next_bp].values[0]])))\n",
    "            # opencv has some issues with dealing with np objects. Cast it manually again\n",
    "            frame = cv2.line(\n",
    "                np.array(frame), coord_0, coord_1, color=(255, 0, 0), thickness=self.line_thickness)\n",
    "            mask = cv2.line(\n",
    "                mask, coord_0, coord_1, color=(255), thickness=self.line_thickness)\n",
    "\n",
    "        # fill out the mask with 1s OUTSIDE of the mask, then invert 0 and 1\n",
    "        # for cv2.floodFill, need a mask that is 2 pixels bigger than the input image\n",
    "        new_mask = np.zeros((mask.shape[0]+2, mask.shape[1]+2), dtype=np.uint8)\n",
    "        cv2.floodFill(mask, new_mask, seedPoint=(0, 0), newVal=124)\n",
    "\n",
    "        final_mask = np.logical_not(new_mask).astype(int)[1:-1, 1:-1]\n",
    "\n",
    "        # ax.imshow(mask)\n",
    "        return {'frame': frame,\n",
    "                'mask': final_mask,\n",
    "                'eyelid_labels_num': len(eyelid_labels)}\n",
    "\n",
    "    def fit_circle_to_pupil(self, frame_num, frame):\n",
    "        \"\"\"\n",
    "        Fit a circle to the pupil if it meets the circle_threshold value\n",
    "        Input:\n",
    "            frame_num: int\n",
    "                A desired frame number\n",
    "            frame: numpy array\n",
    "                A frame to be fitted 3D\n",
    "        Output: dictionary\n",
    "            A dictionary with the fitted frame, center and radius of the fitted circle. If fitting did\n",
    "            not occur, return the original frame with center and raidus as None.\n",
    "            For each key in dictionary:\n",
    "                frame: numpy array\n",
    "                    a numpy array of the frame with pupil circle\n",
    "                center: tuple\n",
    "                    coordinates of the center of the fitted circle. In tuple format\n",
    "                radius: float\n",
    "                    radius of the fitted circle in int format\n",
    "                mask: numpy array\n",
    "                    a binary mask for the fitted circle area\n",
    "                pupil_labels_num: int\n",
    "                    number of pupil labels used for fitting\n",
    "        \"\"\"\n",
    "\n",
    "        mask = np.zeros(frame.shape, dtype=np.uint8)\n",
    "\n",
    "        _, df_x_coords, df_y_coords = self.coords_pcutoff(frame_num)\n",
    "\n",
    "        pupil_labels = [label for label in list(\n",
    "            df_x_coords.index.get_level_values(0)) if 'pupil' in label]\n",
    "\n",
    "        if len(pupil_labels) < self.circle_threshold_num:\n",
    "            # print('Frame number: {} has only 2 or less pupil label. Skip fitting!'.format(\n",
    "            #     frame_num))\n",
    "            center = None\n",
    "            radius = None\n",
    "            final_mask = mask[:, :, 0]\n",
    "\n",
    "        else:\n",
    "            pupil_x = df_x_coords.loc[pupil_labels].values\n",
    "            pupil_y = df_y_coords.loc[pupil_labels].values\n",
    "\n",
    "            pupil_coords = list(zip(pupil_x, pupil_y))\n",
    "\n",
    "            x, y, radius = smallest_enclosing_circle_naive(pupil_coords)\n",
    "\n",
    "            center = (x, y)\n",
    "\n",
    "            # opencv has some issues with dealing with np objects. Cast it manually again\n",
    "            frame = cv2.circle(img=np.array(frame), center=(int(round(x)), int(round(y))),\n",
    "                               radius=int(round(radius)), color=self.circle_color, thickness=self.line_thickness)\n",
    "\n",
    "            mask = cv2.circle(img=mask, center=(int(round(x)), int(round(y))),\n",
    "                              radius=int(round(radius)), color=self.circle_color, thickness=self.line_thickness)\n",
    "\n",
    "            # fill out the mask with 1s OUTSIDE of the mask, then invert 0 and 1\n",
    "            # for cv2.floodFill, need a mask that is 2 pixels bigger than the input image\n",
    "            new_mask = np.zeros(\n",
    "                (mask.shape[0]+2, mask.shape[1]+2), dtype=np.uint8)\n",
    "            cv2.floodFill(mask, new_mask, seedPoint=(0, 0), newVal=1)\n",
    "            final_mask = np.logical_not(new_mask).astype(int)[1:-1, 1:-1]\n",
    "\n",
    "        return {'frame': frame,\n",
    "                'center': center,\n",
    "                'radius': radius,\n",
    "                'mask': final_mask,\n",
    "                'pupil_labels_num': len(pupil_labels)}\n",
    "\n",
    "    def fit_ellipse_to_pupil(self, frame_num, frame):\n",
    "        \"\"\"\n",
    "        Fit an ellipse to pupil iff there exist more than 6 labels. If less than 6, return None\n",
    "        Input:\n",
    "            frame_num: int\n",
    "                A desired frame number\n",
    "            frame: numpy array\n",
    "                A frame to be fitted in 3D\n",
    "        Output: dictionary\n",
    "            A dictionary with the fitted frame, center and radius of the fitted ellipse. If fitting did\n",
    "            not occur, return None.\n",
    "            For each key in dictionary:\n",
    "                frame: numpy array\n",
    "                    a numpy array of the frame with pupil ellipse\n",
    "                center: tuple \n",
    "                    coordinates of the center of the fitted ellipse in tuple of floats\n",
    "                mask: numpy array\n",
    "                    a binary mask for the fitted ellipse area\n",
    "                major_radius: float\n",
    "                    major radius of the fitted ellipse\n",
    "                minor_radius: float\n",
    "                    minor radius of the fitted ellipse\n",
    "                rotation_angle: float\n",
    "                    angle from degree 0 to major_radius                \n",
    "                pupil_labels_num: int\n",
    "                    number of pupil labels used for fitting\n",
    "        \"\"\"\n",
    "\n",
    "        mask = np.zeros(frame.shape, dtype=np.uint8)\n",
    "\n",
    "        _, df_x_coords, df_y_coords = self.coords_pcutoff(frame_num)\n",
    "\n",
    "        pupil_labels = [label for label in list(\n",
    "            df_x_coords.index.get_level_values(0)) if 'pupil' in label]\n",
    "\n",
    "        if len(pupil_labels) < self.ellipse_threshold_num:\n",
    "            # print('Frame number: {} has only 2 or less pupil label. Skip fitting!'.format(\n",
    "            #     frame_num))\n",
    "            center = None\n",
    "            major_radius = None\n",
    "            minor_radius = None\n",
    "            rotation_angle = None\n",
    "            final_mask = mask[:, :, 0]\n",
    "\n",
    "        else:\n",
    "            pupil_x = df_x_coords.loc[pupil_labels].values.round().astype(\n",
    "                np.int32)\n",
    "            pupil_y = df_y_coords.loc[pupil_labels].values.round().astype(\n",
    "                np.int32)\n",
    "\n",
    "            pupil_coords = np.array(\n",
    "                list(zip(pupil_x, pupil_y))).reshape((-1, 1, 2))\n",
    "\n",
    "            # https://docs.opencv.org/2.4/modules/imgproc/doc/structural_analysis_and_shape_descriptors.html#fitellipse\n",
    "            # Python: cv.FitEllipse2(points) → Box2D\n",
    "            # https://docs.opencv.org/3.4.5/db/dd6/classcv_1_1RotatedRect.html\n",
    "            rotated_rect = cv2.fitEllipse(pupil_coords)\n",
    "\n",
    "            # https://docs.opencv.org/3.0-beta/modules/imgproc/doc/drawing_functions.html#ellipse\n",
    "            # cv2.ellipse(img, box, color[, thickness[, lineType]]) → img\n",
    "            frame = cv2.ellipse(np.array(\n",
    "                frame), rotated_rect, color=self.ellipse_color, thickness=self.line_thickness)\n",
    "            mask = cv2.ellipse(np.array(\n",
    "                mask), rotated_rect, color=self.ellipse_color, thickness=self.line_thickness)\n",
    "\n",
    "            # fill out the mask with 1s OUTSIDE of the mask, then invert 0 and 1\n",
    "            # for cv2.floodFill, need a mask that is 2 pixels bigger than the input image\n",
    "            new_mask = np.zeros(\n",
    "                (mask.shape[0]+2, mask.shape[1]+2), dtype=np.uint8)\n",
    "            cv2.floodFill(mask, new_mask, seedPoint=(0, 0), newVal=1)\n",
    "            final_mask = np.logical_not(new_mask).astype(int)[1:-1, 1:-1]\n",
    "\n",
    "            center = rotated_rect[0]\n",
    "            major_radius = rotated_rect[1][1]/2.0\n",
    "            minor_radius = rotated_rect[1][0]/2.0\n",
    "            rotation_angle = rotated_rect[2]\n",
    "\n",
    "        return {'frame': frame,\n",
    "                'center': center,\n",
    "                'mask': final_mask,\n",
    "                'major_radius': major_radius,\n",
    "                'minor_radius': minor_radius,\n",
    "                'rotation_angle': rotation_angle,\n",
    "                'pupil_labels_num': len(pupil_labels)}\n",
    "\n",
    "    def detect_visible_pupil_area(self, eyelid_connect_dict, fit_dict, fitting_method=None):\n",
    "        \"\"\"\n",
    "        Given a frame, find a visible part of the pupil by finding the intersection of pupil and eyelid masks\n",
    "        If pupil mask does not exist(i.e. label < 3), return None\n",
    "        Input:\n",
    "            eyelid_connect_dict: dicionary\n",
    "                An output dictionary from method 'connect_eyelids'\n",
    "            fit_dict: dictionary\n",
    "                An output dictionary from either of the method 'fit_circle_to_pupil' or 'fit_ellipse_to_pupil'\n",
    "            fitting_method: string\n",
    "                A string indicates whether the fitted method was an 'ellipse' or a 'circle'. \n",
    "        Output:\n",
    "            A dictionary that contains the following:\n",
    "                mask: numpy array\n",
    "                    A 3D mask that depicts visible area of pupil. \n",
    "                    If no visible area provided, then it is an np.zeros\n",
    "                visible_portion: signed float\n",
    "                    if visible area exists, then value ranges from 0.0 to 1.0\n",
    "                    if equal to -1.0, not all eyelid labels exist, hence cannot find the visible area\n",
    "                    if equal to -2.0, number of pupil labels do not meet the threshold of the fitting method, hence no fitting performed\n",
    "                    if equal to -3.0, not all eyelid labels exist AND not enough pupil labels to meet the threshold of the fitting method\n",
    "        \"\"\"\n",
    "\n",
    "        color_mask = np.zeros(\n",
    "            shape=[*eyelid_connect_dict['mask'].shape, 3], dtype=np.uint8)\n",
    "\n",
    "        if fitting_method == 'circle':\n",
    "            threshold = self.circle_threshold_num\n",
    "\n",
    "        elif fitting_method == 'ellipse':\n",
    "            threshold = self.ellipse_threshold_num\n",
    "\n",
    "        else:\n",
    "            raise ValueError(\n",
    "                'fitting_method must be provided! It must be either a \"circle\" or an \"ellipse\"')\n",
    "\n",
    "        if fit_dict['pupil_labels_num'] >= threshold and eyelid_connect_dict['eyelid_labels_num'] == 8:\n",
    "\n",
    "            visible_mask = np.logical_and(\n",
    "                fit_dict['mask'], eyelid_connect_dict['mask']).astype(int)\n",
    "\n",
    "            # 126,0,255 for the color\n",
    "            color_mask[visible_mask == 1, 0] = 126\n",
    "            color_mask[visible_mask == 1, 2] = 255\n",
    "\n",
    "            visible_portion = visible_mask.sum() / fit_dict['mask'].sum()\n",
    "\n",
    "        elif fit_dict['pupil_labels_num'] >= threshold and eyelid_connect_dict['eyelid_labels_num'] != 8:\n",
    "            visible_portion = -1.0\n",
    "\n",
    "        elif fit_dict['pupil_labels_num'] < threshold and eyelid_connect_dict['eyelid_labels_num'] == 8:\n",
    "            visible_portion = -2.0\n",
    "\n",
    "        elif fit_dict['pupil_labels_num'] < threshold and eyelid_connect_dict['eyelid_labels_num'] != 8:\n",
    "            visible_portion = -3.0\n",
    "\n",
    "        return dict(mask=color_mask, visible_portion=visible_portion)\n",
    "\n",
    "    def fitted_core(self, frame_num):\n",
    "        \"\"\"\n",
    "        Input:\n",
    "            fig: figure object\n",
    "                Figure object created by configure_plot method\n",
    "            ax: axis object\n",
    "                Axis object created by configure_plot method\n",
    "            frame_num: int\n",
    "                frame number of the video to be analyzed\n",
    "            fitting_method: string\n",
    "                It must be either 'circle' or 'ellipse'. If not provided, raising an error\n",
    "        Output:\n",
    "            A dictionary\n",
    "                'ax_frame': A fitted axis object\n",
    "                'ax_scatter': A scatter axis object that shows where labels are\n",
    "                'ax_mask':\n",
    "        \"\"\"\n",
    "        # it's given in 3 channels but every channel is the same i.e. grayscale\n",
    "\n",
    "        image = self.clip._read_specific_frame(frame_num)\n",
    "\n",
    "        if self.cropped:\n",
    "\n",
    "            x1 = self.cropped_coords[0]\n",
    "            x2 = self.cropped_coords[1]\n",
    "            y1 = self.cropped_coords[2]\n",
    "            y2 = self.cropped_coords[3]\n",
    "\n",
    "            image = image[y1:y2, x1:x2]\n",
    "\n",
    "        eyelid_connected = self.connect_eyelids(\n",
    "            frame_num=frame_num, frame=image)\n",
    "\n",
    "        circle_fit = self.fit_circle_to_pupil(\n",
    "            frame_num=frame_num, frame=eyelid_connected['frame'])\n",
    "\n",
    "        ellipse_fit = self.fit_ellipse_to_pupil(\n",
    "            frame_num=frame_num, frame=eyelid_connected['frame'])\n",
    "\n",
    "        circle_visible = self.detect_visible_pupil_area(eyelid_connect_dict=eyelid_connected,\n",
    "                                                        fit_dict=circle_fit,\n",
    "                                                        fitting_method='circle')\n",
    "\n",
    "        ellipse_visible = self.detect_visible_pupil_area(eyelid_connect_dict=eyelid_connected,\n",
    "                                                         fit_dict=ellipse_fit,\n",
    "                                                         fitting_method='ellipse')\n",
    "\n",
    "        return {'circle_fit': circle_fit,\n",
    "                'ellipse_fit': ellipse_fit,\n",
    "                'circle_visible': circle_visible,\n",
    "                'ellipse_visible': ellipse_visible}\n",
    "\n",
    "    def plot_fitted_frame(self, frame_num, ax=None, fitting_method='circle', save_fig=False):\n",
    "\n",
    "        if ax is None:\n",
    "            fig, ax = self.configure_plot()\n",
    "\n",
    "        # plot bodyparts above the pcutoff\n",
    "        bpindex, x_coords, y_coords = self.coords_pcutoff(frame_num)\n",
    "        ax_scatter = ax.scatter(x_coords.values, y_coords.values, s=self.dotsize**2,\n",
    "                                color=self._label_colors(bpindex), alpha=self.alphavalue)\n",
    "\n",
    "        fitted_core_dict = self.fitted_core(frame_num)\n",
    "\n",
    "        if fitting_method == 'circle':\n",
    "\n",
    "            circle_frame = ax.imshow(fitted_core_dict['circle_fit']['frame'])\n",
    "            circle_mask = ax.imshow(\n",
    "                fitted_core_dict['circle_visible']['mask'], alpha=0.2)\n",
    "\n",
    "        elif fitting_method == 'ellipse':\n",
    "\n",
    "            ellipse_frame = ax.imshow(fitted_core_dict['ellipse_fit']['frame'])\n",
    "            ellipse_mask = ax.imshow(\n",
    "                fitted_core_dict['ellipse_visible']['mask'], alpha=0.2)\n",
    "\n",
    "        else:\n",
    "            raise ValueError(\n",
    "                'fitting method must be either an ellipse or a circle!')\n",
    "\n",
    "        ax.set_title('frame num: ' + str(frame_num), fontsize=self.fontsize)\n",
    "\n",
    "        ax.axis('off')\n",
    "        plt.tight_layout()\n",
    "\n",
    "        if save_fig:\n",
    "            plt.savefig(os.path.join(\n",
    "                self.compressed_cropped_dir_path, 'fitted_frame_' + str(frame_num) + '.png'))\n",
    "\n",
    "        if ax is not None:\n",
    "            return ax\n",
    "\n",
    "    def plot_fitted_multi_frames(self, start, end, fitting_method='circle', save_gif=False):\n",
    "\n",
    "        fig, ax = self.configure_plot()\n",
    "\n",
    "        plt_list = []\n",
    "\n",
    "        for frame_num in range(start, end):\n",
    "\n",
    "            self.plot_fitted_frame(frame_num=frame_num,\n",
    "                                   ax=ax, fitting_method=fitting_method)\n",
    "\n",
    "            fig.canvas.draw()\n",
    "\n",
    "            data = np.fromstring(fig.canvas.tostring_rgb(),\n",
    "                                 dtype=np.uint8, sep='')\n",
    "            data = data.reshape(fig.canvas.get_width_height()[::-1] + (3,))\n",
    "            plt_list.append(data)\n",
    "\n",
    "            display.clear_output(wait=True)\n",
    "            display.display(pl.gcf())\n",
    "            time.sleep(0.5)\n",
    "\n",
    "            plt.cla()\n",
    "\n",
    "        if save_gif:\n",
    "            gif_path = self.video_path.split('.')[0] + '_fitted_' + \\\n",
    "                str(start) + '_' + str(end) + '.gif'\n",
    "            imageio.mimsave(gif_path, plt_list, fps=1)\n",
    "\n",
    "        plt.close('all')\n",
    "\n",
    "    def make_movie(self, start, end, fitting_method='circle'):\n",
    "\n",
    "        import matplotlib.animation as animation\n",
    "\n",
    "        # initlize with start frame\n",
    "        fig, ax = self.configure_plot()\n",
    "\n",
    "        self.plot_fitted_frame(frame_num=start, ax=ax,\n",
    "                               fitting_method=fitting_method)\n",
    "\n",
    "        def update_frame(frame_num):\n",
    "\n",
    "            # clear out the axis\n",
    "            plt.cla()\n",
    "\n",
    "            self.plot_fitted_frame(frame_num=frame_num,\n",
    "                                   ax=ax, fitting_method=fitting_method)\n",
    "\n",
    "        ani = animation.FuncAnimation(fig, update_frame, range(\n",
    "            start+1, end))  # , interval=int(1/self.clip.FPS)\n",
    "        # ani = animation.FuncAnimation(fig, self.plot_fitted_frame, 10)\n",
    "        writer = animation.writers['ffmpeg'](fps=self.clip.FPS)\n",
    "\n",
    "        if self.cropped:\n",
    "            crop_flag = 'cropped'\n",
    "        else:\n",
    "            crop_flag = 'orig'\n",
    "        save_video_path = os.path.join(self.base_dir,\n",
    "                                       '{}_{}_{}_labeled.avi'.format(crop_flag, start, end))\n",
    "\n",
    "        ani.save(save_video_path, writer=writer, dpi=self.dpi)\n",
    "\n",
    "        return ani\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20d9a8ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "class DeeplabcutPlotBodyparts():\n",
    "\n",
    "    def __init__(self, config, bodyparts='all', cropped=False):\n",
    "        \"\"\"\n",
    "        Input:\n",
    "            config: dictionary\n",
    "                A dictionary that contains animal id, session, scan idx, and a path to config\n",
    "            bodyparts: list\n",
    "                A list that contains bodyparts to plot. Each bodypart is in a string format. If none provided,\n",
    "                then by default it plots ALL existing bodyplots in config.yaml file.\n",
    "            cropped: boolean\n",
    "                whether to crop the video or not. Default False\n",
    "            filtering (dict):\n",
    "        \"\"\"\n",
    "\n",
    "        self.config = config\n",
    "\n",
    "        if isinstance(bodyparts, list):\n",
    "            self.bodyparts = bodyparts\n",
    "        else:\n",
    "            self.bodyparts = self.config['bodyparts']\n",
    "\n",
    "        self.cropped = cropped\n",
    "        self.cropped_coords = config['cropped_coords']\n",
    "\n",
    "        self.shuffle = self.config['shuffle']\n",
    "        self.trainingsetindex = self.config['trainingsetindex']\n",
    "\n",
    "        self.project_path = self.config['project_path']\n",
    "        self.orig_video_path = self.config['orig_video_path']\n",
    "        self.base_dir = os.path.dirname(self.orig_video_path)\n",
    "        self.compressed_cropped_dir_path = os.path.join(\n",
    "            os.path.dirname(self.orig_video_path), 'compressed_cropped')\n",
    "        self._case = os.path.basename(config['orig_video_path']).split('.')[0]\n",
    "        self.clip = video_processor.VideoProcessorCV(\n",
    "            fname=self.orig_video_path)\n",
    "\n",
    "        self._trainFraction = self.config['TrainingFraction'][self.trainingsetindex]\n",
    "        self._DLCscorer = auxiliaryfunctions.GetScorerName(\n",
    "            self.config, self.shuffle, self._trainFraction)\n",
    "\n",
    "        self.label_path = os.path.join(\n",
    "            self.compressed_cropped_dir_path, self._case + '_compressed_cropped' + self._DLCscorer + '.h5')\n",
    "\n",
    "        self.df_label = pd.read_hdf(self.label_path)\n",
    "\n",
    "        self.df_bodyparts = self.df_label[self._DLCscorer][self.bodyparts]\n",
    "\n",
    "        self.df_bodyparts_likelihood = self.df_bodyparts.iloc[:, self.df_bodyparts.columns.get_level_values(\n",
    "            1) == 'likelihood']\n",
    "\n",
    "        self.df_bodyparts_x = self.df_bodyparts.iloc[:,\n",
    "                                                     self.df_bodyparts.columns.get_level_values(1) == 'x']\n",
    "        self.df_bodyparts_y = self.df_bodyparts.iloc[:,\n",
    "                                                     self.df_bodyparts.columns.get_level_values(1) == 'y']\n",
    "\n",
    "        # in mm. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3310398/#R13\n",
    "        self._pupil_diameter = 3.0 \n",
    "\n",
    "        # obtain median left to right eyelid distance\n",
    "        right = self.df_bodyparts['eyelid_right'].values[:, :2]\n",
    "        left = self.df_bodyparts['eyelid_left'].values[:, :2]\n",
    "\n",
    "        self.median_left_right = np.median(\n",
    "            np.sqrt(np.einsum('ij,ij->i', left-right, left-right)))\n",
    "\n",
    "        # obtain pixel to diameter ratio\n",
    "        self._pixel_to_diameter_ratio = self.median_left_right/self._pupil_diameter\n",
    "\n",
    "        if not self.cropped:\n",
    "\n",
    "            self.nx = self.clip.width()\n",
    "            self.ny = self.clip.height()\n",
    "            self.df_bodyparts_x = self.df_bodyparts.iloc[:,\n",
    "                                                         self.df_bodyparts.columns.get_level_values(1) == 'x'] + self.cropped_coords[0]\n",
    "            self.df_bodyparts_y = self.df_bodyparts.iloc[:,\n",
    "                                                         self.df_bodyparts.columns.get_level_values(1) == 'y'] + self.cropped_coords[2]\n",
    "\n",
    "        else:\n",
    "\n",
    "            if self.cropped_coords is None:\n",
    "                raise ValueError(\n",
    "                    \"cropped_coords are not provided! Must be in list with 4 elmnts long!\")\n",
    "\n",
    "            if len(self.cropped_coords) != 4:\n",
    "                raise ValueError(\n",
    "                    \"Only provided {} coordinates! U need 4!\".format(len(self.cropped_coords)))\n",
    "\n",
    "            # self.df_bodyparts_x = self.df_bodyparts.iloc[:,\n",
    "            #                                             self.df_bodyparts.columns.get_level_values(1) == 'x'] - self.cropped_coords[0]\n",
    "            # self.df_bodyparts_y = self.df_bodyparts.iloc[:,\n",
    "            #                                             self.df_bodyparts.columns.get_level_values(1) == 'y'] - self.cropped_coords[2]\n",
    "\n",
    "            self.nx = self.clip.width() - self.cropped_coords[0]\n",
    "            self.ny = self.clip.height() - self.cropped_coords[2]\n",
    "\n",
    "        # plotting properties\n",
    "        self._dotsize = 7\n",
    "        self._line_thickness = 1\n",
    "        self._pcutoff = self.config['pcutoff']\n",
    "        self._colormap = self.config['colormap']\n",
    "        self._label_colors = plotting.get_cmap(\n",
    "            len(self.bodyparts), name=self._colormap)\n",
    "        self._alphavalue = self.config['alphavalue']\n",
    "        self._fig_size = [12, 8]\n",
    "        self._dpi = 100\n",
    "        self._fontsize = 30\n",
    "\n",
    "        self.tf_likelihood_array = self.df_bodyparts_likelihood.values > self._pcutoff\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10a96ef2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def smallest_enclosing_circle_naive(points):\n",
    "    # Degenerate cases\n",
    "    if len(points) == 0:\n",
    "        return None\n",
    "    elif len(points) == 1:\n",
    "        return (points[0][0], points[0][1], 0)\n",
    "\n",
    "    result = None\n",
    "    # Try all unique triples\n",
    "    for i in range(len(points)):\n",
    "        p = points[i]\n",
    "        for j in range(i + 1, len(points)):\n",
    "            q = points[j]\n",
    "            for k in range(j + 1, len(points)):\n",
    "                r = points[k]\n",
    "                c = make_circumcircle(p, q, r)\n",
    "                if c is not None and (result is None or c[2] < result[2]) and \\\n",
    "                        all(is_in_circle(c, s) for s in points):\n",
    "                    result = c\n",
    "\n",
    "    if result is None:\n",
    "        raise AssertionError()\n",
    "    return result\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90f41e0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def make_circumcircle(a, b, c):\n",
    "    # Mathematical algorithm from Wikipedia: Circumscribed circle\n",
    "    ox = (min(a[0], b[0], c[0]) + max(a[0], b[0], c[0])) / 2.0\n",
    "    oy = (min(a[1], b[1], c[1]) + max(a[1], b[1], c[1])) / 2.0\n",
    "    ax = a[0] - ox\n",
    "    ay = a[1] - oy\n",
    "    bx = b[0] - ox\n",
    "    by = b[1] - oy\n",
    "    cx = c[0] - ox\n",
    "    cy = c[1] - oy\n",
    "    d = (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by)) * 2.0\n",
    "    if d == 0.0:\n",
    "        return None\n",
    "    x = ox + ((ax*ax + ay*ay) * (by - cy) + (bx*bx + by*by)\n",
    "              * (cy - ay) + (cx*cx + cy*cy) * (ay - by)) / d\n",
    "    y = oy + ((ax*ax + ay*ay) * (cx - bx) + (bx*bx + by*by)\n",
    "              * (ax - cx) + (cx*cx + cy*cy) * (bx - ax)) / d\n",
    "    ra = math.hypot(x - a[0], y - a[1])\n",
    "    rb = math.hypot(x - b[0], y - b[1])\n",
    "    rc = math.hypot(x - c[0], y - c[1])\n",
    "    return (x, y, max(ra, rb, rc))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72bfa097",
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_in_circle(c, p):\n",
    "    return c is not None and math.hypot(p[0] - c[0], p[1] - c[1]) <= c[2] * _MULTIPLICATIVE_EPSILON\n"
   ]
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}