[dbf742]: / Inference_Video.ipynb

Download this file

656 lines (655 with data), 28.9 kB

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import necessary packages (run once upon startup)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using TensorFlow backend.\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division \n",
    "import os\n",
    "import pandas as pd\n",
    "from pandas import ExcelWriter\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use(\"ggplot\")\n",
    "%matplotlib inline\n",
    "\n",
    "\n",
    "from skimage.transform import resize\n",
    "from skimage.morphology import skeletonize\n",
    "from scipy.signal import resample, savgol_filter, butter, filtfilt\n",
    "from PIL import Image, ImageDraw\n",
    "import cv2\n",
    "\n",
    "import tensorflow as tf\n",
    "\n",
    "from keras import backend as K\n",
    "from keras.models import Model, load_model\n",
    "from keras.preprocessing.image import ImageDataGenerator, array_to_img, img_to_array, load_img"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define custom functions and load trained models (run once upon startup)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Intersection over union (IoU), a measure of labelling accuracy (sometimes also called Jaccard score)\n",
    "def IoU(y_true, y_pred, smooth=1):\n",
    "    intersection = K.sum(K.abs(y_true * y_pred), axis=-1)\n",
    "    union = K.sum(y_true,-1) + K.sum(y_pred,-1) - intersection\n",
    "    iou = (intersection + smooth) / ( union + smooth)\n",
    "    return iou\n",
    "\n",
    "# Function to sort contours from proximal to distal (the bounding boxes are not used)\n",
    "def sort_contours(cnts):\n",
    "    # initialize the reverse flag and sort index\n",
    "    i = 1\n",
    "    # construct the list of bounding boxes and sort them from top to bottom\n",
    "    boundingBoxes = [cv2.boundingRect(c) for c in cnts]\n",
    "    (cnts, boundingBoxes) = zip(*sorted(zip(cnts, boundingBoxes), key=lambda b:b[1][i], reverse=False))\n",
    " \n",
    "    return (cnts, boundingBoxes)\n",
    "\n",
    "# Find only the coordinates representing one edge of a contour. edge: T (top) or B (bottom)\n",
    "def contour_edge(edge, contour):\n",
    "    pts = list(contour)\n",
    "    ptsT = sorted(pts, key=lambda k: [k[0][0], k[0][1]])\n",
    "    allx = []\n",
    "    ally = []\n",
    "    for a in range(0,len(ptsT)):\n",
    "        allx.append(ptsT[a][0,0])\n",
    "        ally.append(ptsT[a][0,1])\n",
    "    un = np.unique(allx)\n",
    "    #sumA = 0\n",
    "    leng = len(un)-1\n",
    "    x = []\n",
    "    y = []\n",
    "    for each in range(5,leng-5): # Ignore 1st and last 5 points to avoid any curves\n",
    "        indices = [i for i, x in enumerate(allx) if x == un[each]]\n",
    "        if edge == 'T':\n",
    "            loc = indices[0]\n",
    "        else:\n",
    "            loc = indices[-1]\n",
    "        x.append(ptsT[loc][0,0])\n",
    "        y.append(ptsT[loc][0,1])\n",
    "    return np.array(x),np.array(y)\n",
    "\n",
    "def intersection(L1, L2):\n",
    "    D  = L1[0] * L2[1] - L1[1] * L2[0]\n",
    "    Dx = L1[2] * L2[1] - L1[1] * L2[2]\n",
    "    Dy = L1[0] * L2[2] - L1[2] * L2[0]\n",
    "    if D != 0:\n",
    "        x = Dx / D\n",
    "        y = Dy / D\n",
    "        return x,y\n",
    "    else:\n",
    "        return False\n",
    "\n",
    "# Function to detect mouse clicks for the purpose of image calibration\n",
    "def mclick(event, x, y, flags, param):\n",
    "    # grab references to the global variables\n",
    "    global mlocs\n",
    "\n",
    "    # if the left mouse button was clicked, record the (x, y) coordinates\n",
    "    if event == cv2.EVENT_LBUTTONDOWN:\n",
    "        mlocs.append(y)\n",
    "        \n",
    "# Function to compute the distance between 2 x,y points\n",
    "def distFunc(x1, y1, x2, y2):\n",
    "    xdist = (x2 - x1)**2\n",
    "    ydist = (y2 - y1)**2\n",
    "    return np.sqrt(xdist + ydist)\n",
    "\n",
    "###############################################################################\n",
    "\n",
    "# IMPORT THE TRAINED MODELS\n",
    "\n",
    "# load the aponeurosis model\n",
    "model_apo = load_model('./models/model-apo2-nc.h5', custom_objects={'IoU': IoU})\n",
    "\n",
    "# load the fascicle model\n",
    "modelF = load_model('./models/model-fascSnippets2-nc.h5', custom_objects={'IoU': IoU})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import Video that you want to analyse (set location under 'vpath')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DEFINE THE PATH OF YOUR VIDEO \n",
    "vpath = 'D:/Unet annotations/UltratrackCompVideos/Olivier_vids/GM_walk/sp1_10.avi'\n",
    "# vpath = 'D:/Unet annotations/UltratrackCompVideos/MG_MVC.avi'\n",
    "\n",
    "# Video properties (do not edit)\n",
    "cap = cv2.VideoCapture(vpath)\n",
    "vid_len = int(cap.get(cv2.CAP_PROP_FRAME_COUNT))\n",
    "vid_width = int(cap.get(cv2.CAP_PROP_FRAME_WIDTH))\n",
    "vid_height = int(cap.get(cv2.CAP_PROP_FRAME_HEIGHT))\n",
    "vid_fps = int(cap.get(cv2.CAP_PROP_FPS))\n",
    "count = 0\n",
    "dataout = []\n",
    "indices = [i for i, a in enumerate(vpath) if a == '/']\n",
    "dots = [i for i, a in enumerate(vpath) if a == '.']\n",
    "filename = './analysedVideos/' + vpath[indices[-1]+1:dots[-1]] + '.avi' # indices\n",
    "vid_out = cv2.VideoWriter(filename,cv2.VideoWriter_fourcc(*'DIVX'), vid_fps, (vid_width, vid_height))\n",
    "calibDist = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define settings\n",
    "apo_threshold = 0.3                     # Sensitivity threshold for detecting aponeuroses\n",
    "fasc_threshold = 0.10                   # Sensitivity threshold for detecting fascicles\n",
    "fasc_cont_thresh = 40                   # Minimum accepted contour length for fascicles (px)\n",
    "flip = 0                                # If fascicles are oriented bottom-left to top-right, leave as 0. Otherwise set to 1\n",
    "min_width = 40                          # Minimum acceptable distance between aponeuroses\n",
    "curvature = 1                           # Set to 3 for curved fascicles or 1 for a straight line\n",
    "min_pennation = 10                      # Minimum and maximum acceptable pennation angles\n",
    "max_pennation = 40"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "list index out of range",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-19-a408670d9994>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     17\u001b[0m         \u001b[0mcv2\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdestroyAllWindows\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     18\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 19\u001b[1;33m     \u001b[0mcalibDist\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmlocs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mmlocs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     20\u001b[0m     \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mstr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mspacing\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m+\u001b[0m \u001b[1;34m' mm corresponds to '\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mstr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcalibDist\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m+\u001b[0m \u001b[1;34m' pixels'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mIndexError\u001b[0m: list index out of range"
     ]
    }
   ],
   "source": [
    "# OPTIONAL\n",
    "# Calibrate the analysis by clicking on 2 points in the image, followed by the 'q' key. These two points should be 1cm apart\n",
    "# Alternatively, change the spacing setting below\n",
    "# NOTE: Here we assume that the points are spaced apart in the y/vertical direction of the image\n",
    "spacing = 10.0 # Space between the two calibration markers (mm)\n",
    "mlocs = []\n",
    "for cal in range(0,1):\n",
    "    _,frame = cap.read()\n",
    "\n",
    "    # display the image and wait for a keypress\n",
    "    cv2.imshow(\"image\", frame)\n",
    "    cv2.setMouseCallback(\"image\", mclick)\n",
    "    key = cv2.waitKey(0)\n",
    " \n",
    "    # if the 'q' key is pressed, break from the loop\n",
    "    if key == ord(\"q\"):\n",
    "        cv2.destroyAllWindows()\n",
    "\n",
    "    calibDist.append(np.abs(mlocs[0] - mlocs[1]))\n",
    "    print(str(spacing) + ' mm corresponds to ' + str(calibDist[0]) + ' pixels')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\neilj\\AppData\\Local\\Continuum\\anaconda3\\envs\\tf_gpu\\lib\\site-packages\\skimage\\transform\\_warps.py:110: UserWarning: Anti-aliasing will be enabled by default in skimage 0.15 to avoid aliasing artifacts when down-sampling images.\n",
      "  warn(\"Anti-aliasing will be enabled by default in skimage 0.15 to \"\n",
      "C:\\Users\\neilj\\AppData\\Local\\Continuum\\anaconda3\\envs\\tf_gpu\\lib\\site-packages\\skimage\\transform\\_warps.py:105: UserWarning: The default mode, 'constant', will be changed to 'reflect' in skimage 0.15.\n",
      "  warn(\"The default mode, 'constant', will be changed to 'reflect' in \"\n"
     ]
    }
   ],
   "source": [
    "# Analyse whole video\n",
    "fasc_l_all = []\n",
    "pennation_all = []\n",
    "x_lows_all = []\n",
    "x_highs_all = []\n",
    "thickness_all = []\n",
    "\n",
    "for a in range(0, vid_len-1):\n",
    "# for a in range(0, 10):\n",
    "\n",
    "    # FORMAT EACH FRAME, RESHAPE AND COMPUTE NN PREDICTIONS\n",
    "    _,frame = cap.read()\n",
    "    img = img_to_array(frame)\n",
    "    if flip == 1:\n",
    "        img = np.fliplr(img)\n",
    "    img_orig = img # Make a copy\n",
    "    img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n",
    "    h = img.shape[0]\n",
    "    w = img.shape[1]\n",
    "    img = np.reshape(img,[-1, h, w,1])\n",
    "    img = resize(img, (1, 512, 512, 1), mode = 'constant', preserve_range = True)\n",
    "    img = img/255.0\n",
    "    \n",
    "    pred_apo = model_apo.predict(img)\n",
    "    pred_apo_t = (pred_apo > apo_threshold).astype(np.uint8) \n",
    "    pred_fasc = modelF.predict(img)\n",
    "    pred_fasc_t = (pred_fasc > fasc_threshold).astype(np.uint8) \n",
    "\n",
    "    img = resize(img, (1, h, w,1))\n",
    "    img = np.reshape(img, (h, w))\n",
    "    pred_apo = resize(pred_apo, (1, h, w,1))\n",
    "    pred_apo = np.reshape(pred_apo, (h, w))\n",
    "    pred_apo_t = resize(pred_apo_t, (1, h, w,1))\n",
    "    pred_apo_t = np.reshape(pred_apo_t, (h, w))\n",
    "    pred_fasc = resize(pred_fasc, (1, h, w,1))\n",
    "    pred_fasc = np.reshape(pred_fasc, (h, w))\n",
    "    pred_fasc_t = resize(pred_fasc_t, (1, h, w,1))\n",
    "    pred_fasc_t = np.reshape(pred_fasc_t, (h, w))\n",
    "    \n",
    "    #############################################\n",
    "    \n",
    "    # COMPUTE CONTOURS TO IDENTIFY THE APONEUROSES\n",
    "    _, thresh = cv2.threshold(pred_apo_t, 0, 255, cv2.THRESH_BINARY) # Or set lower threshold above and use pred_ind_t here\n",
    "    thresh = thresh.astype('uint8')\n",
    "    contours, hierarchy = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)\n",
    "    \n",
    "#     contours_re = []\n",
    "#     for contour in contours: # Remove any contours that are very small\n",
    "#         if len(contour) > 600:\n",
    "#             contours_re.append(contour)\n",
    "    contours = [i for i in contours if len(i) > 600]\n",
    "#     contours = contours_re\n",
    "    contours,_ = sort_contours(contours) # Sort contours from top to bottom\n",
    "    \n",
    "#     mask_apo = np.zeros(thresh.shape,np.uint8)\n",
    "    contours_re2 = []\n",
    "    for contour in contours:\n",
    "#         cv2.drawContours(mask_apo,[contour],0,255,-1)\n",
    "        pts = list(contour)\n",
    "        ptsT = sorted(pts, key=lambda k: [k[0][0], k[0][1]]) # Sort each contour based on x values\n",
    "        allx = []\n",
    "        ally = []\n",
    "        for aa in range(0,len(ptsT)):\n",
    "            allx.append(ptsT[aa][0,0])\n",
    "            ally.append(ptsT[aa][0,1])\n",
    "        app = np.array(list(zip(allx,ally)))\n",
    "        contours_re2.append(app)\n",
    "\n",
    "    # Merge nearby contours\n",
    "#     countU = 0\n",
    "    xs1 = []\n",
    "    xs2 = []\n",
    "    ys1 = []\n",
    "    ys2 = []\n",
    "    maskT = np.zeros(thresh.shape,np.uint8)\n",
    "    for cnt in contours_re2:\n",
    "        ys1.append(cnt[0][1])\n",
    "        ys2.append(cnt[-1][1])\n",
    "        xs1.append(cnt[0][0])\n",
    "        xs2.append(cnt[-1][0])\n",
    "        cv2.drawContours(maskT,[cnt],0,255,-1)\n",
    "\n",
    "    for countU in range(0,len(contours_re2)-1):\n",
    "        if xs1[countU+1] > xs2[countU]: # Check if x of contour2 is higher than x of contour 1\n",
    "            y1 = ys2[countU]\n",
    "            y2 = ys1[countU+1]\n",
    "            if y1-10 <= y2 <= y1+10:\n",
    "                m = np.vstack((contours_re2[countU], contours_re2[countU+1]))\n",
    "                cv2.drawContours(maskT,[m],0,255,-1)\n",
    "        countU += 1\n",
    "\n",
    "    maskT[maskT > 0] = 1\n",
    "    skeleton = skeletonize(maskT).astype(np.uint8)\n",
    "    kernel = np.ones((3,7), np.uint8) \n",
    "    dilate = cv2.dilate(skeleton, kernel, iterations=15)\n",
    "    erode = cv2.erode(dilate, kernel, iterations=10)\n",
    "\n",
    "    contoursE, hierarchy = cv2.findContours(erode, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)\n",
    "    mask_apoE = np.zeros(thresh.shape,np.uint8)\n",
    "#     contours_reE = []\n",
    "#     for contour in contoursE: # Remove any contours that are very small\n",
    "#         if len(contour) > 600:\n",
    "#             contours_reE.append(contour)\n",
    "#     contoursE = contours_reE\n",
    "    contoursE = [i for i in contoursE if len(i) > 600]\n",
    "    for contour in contoursE:\n",
    "        cv2.drawContours(mask_apoE,[contour],0,255,-1)\n",
    "    contoursE,_ = sort_contours(contoursE)\n",
    "\n",
    "    if len(contoursE) >= 2:\n",
    "    \n",
    "        # Get the x,y coordinates of the upper/lower edge of the 2 aponeuroses\n",
    "        # NOTE: THE APONEUROSES DISPLAYED IN THE IMAGES ARE FROM mask_apoE, NOT THE CONTOUR EDGES BELOW\n",
    "        upp_x,upp_y = contour_edge('B', contoursE[0])\n",
    "        if contoursE[1][0,0,1] > (contoursE[0][0,0,1] + min_width):\n",
    "            low_x,low_y = contour_edge('T', contoursE[1])\n",
    "        else:\n",
    "            low_x,low_y = contour_edge('T', contoursE[2])\n",
    "\n",
    "        upp_y_new = savgol_filter(upp_y, 81, 2) # window size, polynomial order\n",
    "        low_y_new = savgol_filter(low_y, 81, 2)\n",
    "\n",
    "        # Make a binary mask to only include fascicles within the region between the 2 aponeuroses\n",
    "        ex_mask = np.zeros(thresh.shape,np.uint8)\n",
    "        ex_1 = 0\n",
    "        ex_2 = np.minimum(len(low_x), len(upp_x))\n",
    "        for ii in range(ex_1, ex_2):\n",
    "            ymin = int(np.floor(upp_y_new[ii]))\n",
    "            ymax = int(np.ceil(low_y_new[ii]))\n",
    "\n",
    "            ex_mask[:ymin, ii] = 0\n",
    "            ex_mask[ymax:, ii] = 0\n",
    "            ex_mask[ymin:ymax, ii] = 255\n",
    "\n",
    "        # Calculate slope of central portion of each aponeurosis & use this to compute muscle thickness\n",
    "        Alist = list(set(upp_x).intersection(low_x))\n",
    "        Alist = sorted(Alist)\n",
    "        Alen = len(list(set(upp_x).intersection(low_x))) # How many values overlap between x-axes\n",
    "        A1 = int(Alist[0] + (.33 * Alen))\n",
    "        A2 = int(Alist[0] + (.66 * Alen)) \n",
    "        mid = int((A2-A1) / 2 + A1)\n",
    "        mindist = 10000\n",
    "        upp_ind = np.where(upp_x==mid)\n",
    "\n",
    "        if upp_ind == len(upp_x):\n",
    "            upp_ind -= 1\n",
    "\n",
    "        for val in range(A1, A2):\n",
    "            if val >= len(low_x):\n",
    "                continue\n",
    "            else:\n",
    "                dist = distFunc(upp_x[upp_ind], upp_y_new[upp_ind], low_x[val], low_y_new[val])\n",
    "                if dist < mindist:\n",
    "                    mindist = dist\n",
    "\n",
    "        # Add aponeuroses to a mask for display\n",
    "        imgT = np.zeros((h,w,3), np.uint8)\n",
    "\n",
    "        # Compute functions to approximate the shape of the aponeuroses\n",
    "        zUA = np.polyfit(upp_x, upp_y_new, 2) # 2nd order polynomial\n",
    "        g = np.poly1d(zUA)\n",
    "        zLA = np.polyfit(low_x, low_y_new, 2)\n",
    "        h = np.poly1d(zLA)\n",
    "\n",
    "        mid = (low_x[-1]-low_x[0])/2 + low_x[0] # Find middle of the aponeurosis\n",
    "        x1 = np.linspace(low_x[0]-700, low_x[-1]+700, 10000) # Extrapolate polynomial fits to either side of the mid-point\n",
    "        y_UA = g(x1)\n",
    "        y_LA = h(x1)\n",
    "\n",
    "        new_X_UA = np.linspace(mid-700, mid+700, 5000) # Extrapolate x,y data using f function\n",
    "        new_Y_UA = g(new_X_UA)\n",
    "        new_X_LA = np.linspace(mid-700, mid+700, 5000) # Extrapolate x,y data using f function\n",
    "        new_Y_LA = h(new_X_LA)\n",
    "\n",
    "        #############################################\n",
    "\n",
    "        # COMPUTE CONTOURS TO IDENTIFY FASCICLES/FASCICLE ORIENTATION\n",
    "        _, threshF = cv2.threshold(pred_fasc_t, 0, 255, cv2.THRESH_BINARY) \n",
    "        threshF = threshF.astype('uint8')\n",
    "        contoursF, hierarchy = cv2.findContours(threshF, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)\n",
    "\n",
    "        # Remove any contours that are very small\n",
    "#         contours_re = []\n",
    "        maskF = np.zeros(threshF.shape,np.uint8)\n",
    "        for contour in contoursF: # Remove any contours that are very small\n",
    "            if len(contour) > fasc_cont_thresh:\n",
    "#                 contours_re.append(contour)\n",
    "                cv2.drawContours(maskF,[contour],0,255,-1)       \n",
    "\n",
    "        # Only include fascicles within the region of the 2 aponeuroses  \n",
    "        mask_Fi = maskF & ex_mask \n",
    "        contoursF2, hierarchy = cv2.findContours(mask_Fi, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)\n",
    "\n",
    "        contoursF3 = []\n",
    "        for contour in contoursF2:\n",
    "            if len(contour) > fasc_cont_thresh:\n",
    "                contoursF3.append(contour)\n",
    "\n",
    "        xs = []\n",
    "        ys = []\n",
    "        fas_ext = []\n",
    "        fasc_l = []\n",
    "        pennation = []\n",
    "        x_low1 = []\n",
    "        x_high1 = []\n",
    "#         counter = 0\n",
    "\n",
    "        for cnt in contoursF3:\n",
    "    #         x,y = contour_edge('B', contoursF2[counter])\n",
    "            x,y = contour_edge('B', cnt)\n",
    "            if len(x) == 0:\n",
    "                continue\n",
    "            z = np.polyfit(np.array(x), np.array(y), 1)\n",
    "            f = np.poly1d(z)\n",
    "            newX = np.linspace(-400, w+400, 5000) # Extrapolate x,y data using f function\n",
    "            newY = f(newX)\n",
    "\n",
    "            # Find intersection between each fascicle and the aponeuroses.\n",
    "            diffU = newY-new_Y_UA # Find intersections\n",
    "            locU = np.where(diffU == min(diffU, key=abs))[0]\n",
    "            diffL = newY-new_Y_LA\n",
    "            locL = np.where(diffL == min(diffL, key=abs))[0]\n",
    "\n",
    "            coordsX = newX[int(locL):int(locU)]\n",
    "            coordsY = newY[int(locL):int(locU)]\n",
    "\n",
    "            if locL >= 4950:\n",
    "                Apoangle = int(np.arctan((new_Y_LA[locL-50]-new_Y_LA[locL-50])/(new_X_LA[locL]-new_X_LA[locL-50]))*180/np.pi)\n",
    "            else:\n",
    "                Apoangle = int(np.arctan((new_Y_LA[locL]-new_Y_LA[locL+50])/(new_X_LA[locL+50]-new_X_LA[locL]))*180/np.pi) # Angle relative to horizontal\n",
    "            Apoangle = 90 + abs(Apoangle)\n",
    "\n",
    "            # Don't include fascicles that are completely outside of the field of view or\n",
    "            # those that don't pass through central 1/3 of the image\n",
    "    #         if np.sum(coordsX) > 0 and coordsX[-1] > 0 and coordsX[0] < np.maximum(upp_x[-1],low_x[-1]) and coordsX[-1] - coordsX[0] < w:\n",
    "            if np.sum(coordsX) > 0 and coordsX[-1] > 0 and coordsX[0] < np.maximum(upp_x[-1],low_x[-1]) and Apoangle != float('nan'):\n",
    "                FascAng = float(np.arctan((coordsX[0]-coordsX[-1])/(new_Y_LA[locL]-new_Y_UA[locU]))*180/np.pi)*-1\n",
    "                ActualAng = Apoangle-FascAng\n",
    "\n",
    "                if ActualAng <= max_pennation and ActualAng >= min_pennation: # Don't include 'fascicles' beyond a range of pennation angles\n",
    "                    length1 = np.sqrt((newX[locU] - newX[locL])**2 + (y_UA[locU] - y_LA[locL])**2)\n",
    "                    fasc_l.append(length1[0]) # Calculate fascicle length\n",
    "                    pennation.append(Apoangle-FascAng)\n",
    "                    x_low1.append(coordsX[0].astype('int32'))\n",
    "                    x_high1.append(coordsX[-1].astype('int32'))\n",
    "                    coords = np.array(list(zip(coordsX.astype('int32'), coordsY.astype('int32'))))\n",
    "                    cv2.polylines(imgT, [coords], False, (20, 15, 200), 3)\n",
    "#             counter += 1 \n",
    "\n",
    "        # Store the results for each frame and normalise using scale factor (if calibration was done above)\n",
    "        try:\n",
    "            midthick = mindist[0] # Muscle thickness\n",
    "        except:\n",
    "            midthick = mindist\n",
    "\n",
    "        if 'calibDist' in locals() and len(calibDist) >0:\n",
    "            fasc_l = fasc_l / (calibDist[0]/10)\n",
    "            midthick = midthick / (calibDist[0]/10)\n",
    "            \n",
    "    else:\n",
    "        fasc_l = []\n",
    "        pennation = []\n",
    "        x_low1 = []\n",
    "        x_high1 = []\n",
    "        imgT = np.zeros((h,w,3), np.uint8)\n",
    "        fasc_l.append(float(\"nan\"))\n",
    "        pennation.append(float(\"nan\"))\n",
    "        x_low1.append(float(\"nan\"))\n",
    "        x_high1.append(float(\"nan\"))\n",
    "        midthick = float(\"nan\")\n",
    "    \n",
    "    fasc_l_all.append(fasc_l)\n",
    "    pennation_all.append(pennation)\n",
    "    x_lows_all.append(x_low1)\n",
    "    x_highs_all.append(x_high1)\n",
    "    \n",
    "    thickness_all.append(midthick)\n",
    "    \n",
    "    #############################################\n",
    "    \n",
    "    # DISPLAY EACH PROCESSED IMAGE AND METRICS\n",
    "    \n",
    "    img_orig[mask_apoE > 0] = (235, 25, 42)\n",
    "    \n",
    "    comb = cv2.addWeighted(img_orig.astype(np.uint8), 1, imgT, 0.8, 0) \n",
    "    vid_out.write(comb) # Write each image to video file\n",
    "    cv2.putText(comb, ('Frame: ' + str(a+1) + ' of ' + str(vid_len)), (125,350), cv2.FONT_HERSHEY_DUPLEX, 1, (249,249,249))\n",
    "    cv2.putText(comb, ('Pennation angle: ' + str('%.1f' % np.median(pennation_all[-1])) + ' deg'), (125,410), cv2.FONT_HERSHEY_DUPLEX, 1, (249,249,249))\n",
    "    if calibDist:\n",
    "        cv2.putText(comb, ('Fascicle length: ' + str('%.2f' % np.median(fasc_l_all[-1]) + ' mm')), (125,380), cv2.FONT_HERSHEY_DUPLEX, 1, (249,249,249))\n",
    "        cv2.putText(comb, ('Thickness at centre: ' + str('%.1f' % thickness_all[-1]) + ' mm'), (125,440), cv2.FONT_HERSHEY_DUPLEX, 1, (249,249,249))\n",
    "    else:\n",
    "        cv2.putText(comb, ('Fascicle length: ' + str('%.2f' % np.median(fasc_l_all[-1]) + ' px')), (125,380), cv2.FONT_HERSHEY_DUPLEX, 1, (249,249,249))\n",
    "        cv2.putText(comb, ('Thickness at centre: ' + str('%.1f' % thickness_all[-1]) + ' px'), (125,440), cv2.FONT_HERSHEY_DUPLEX, 1, (249,249,249))\n",
    "    \n",
    "    cv2.imshow('Analysed image', comb)\n",
    "    \n",
    "#     count += 1       \n",
    "\n",
    "    if cv2.waitKey(10) & 0xFF == ord('q'): # Press 'q' to stop the analysis\n",
    "        break\n",
    "\n",
    "cap.release()\n",
    "vid_out.release()\n",
    "cv2.destroyAllWindows() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "## Prepare data and export to excel file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "fl = np.zeros([len(fasc_l_all),len(max(fasc_l_all, key = lambda x: len(x)))])\n",
    "pe = np.zeros([len(pennation_all),len(max(pennation_all, key = lambda x: len(x)))])\n",
    "xl = np.zeros([len(x_lows_all),len(max(x_lows_all, key = lambda x: len(x)))])\n",
    "xh = np.zeros([len(x_highs_all),len(max(x_highs_all, key = lambda x: len(x)))])\n",
    "\n",
    "for i,j in enumerate(fasc_l_all):\n",
    "    fl[i][0:len(j)] = j\n",
    "fl[fl==0] = np.nan\n",
    "for i,j in enumerate(pennation_all):\n",
    "    pe[i][0:len(j)] = j\n",
    "pe[pe==0] = np.nan\n",
    "for i,j in enumerate(x_lows_all):\n",
    "    xl[i][0:len(j)] = j\n",
    "xl[xl==0] = np.nan\n",
    "for i,j in enumerate(x_highs_all):\n",
    "    xh[i][0:len(j)] = j\n",
    "xh[xh==0] = np.nan\n",
    "\n",
    "df1 = pd.DataFrame(data=fl)\n",
    "df2 = pd.DataFrame(data=pe)\n",
    "df3 = pd.DataFrame(data=xl)\n",
    "df4 = pd.DataFrame(data=xh)\n",
    "df5 = pd.DataFrame(data=thickness_all)\n",
    "\n",
    "# Create a Pandas Excel writer and your xlsx filename (same as the video filename)\n",
    "# writer = ExcelWriter('./analysedVideos/' + os.path.splitext(filename)[0] + '.xlsx')\n",
    "writer = ExcelWriter('./analysedVideos/' + vpath[indices[-1]+1:dots[-1]] + '.xlsx')\n",
    "\n",
    "# Write each dataframe to a different worksheet.\n",
    "df1.to_excel(writer, sheet_name='Fasc_length')\n",
    "df2.to_excel(writer, sheet_name='Pennation')\n",
    "df3.to_excel(writer, sheet_name='X_low')\n",
    "df4.to_excel(writer, sheet_name='X_high')\n",
    "df5.to_excel(writer, sheet_name='Thickness')\n",
    "\n",
    "# Close the Pandas Excel writer and output the Excel file\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}