kymo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "96d5f190-f8ec-43bc-b149-e9b8152d0de2",
   "metadata": {},
   "outputs": [],
   "source": [
    "##Import libraries\n",
    "import torch\n",
    "from torch import nn\n",
    "import os\n",
    "import h5py\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "from matplotlib import cm\n",
    "import matplotlib.image as mpimg\n",
    "import numpy as np\n",
    "import cv2\n",
    "import math\n",
    "import time\n",
    "from scipy.interpolate import UnivariateSpline, interp1d\n",
    "from scipy import interpolate\n",
    "from scipy import optimize\n",
    "from scipy import integrate\n",
    "from scipy import signal\n",
    "#from arc_length_1 import *\n",
    "from scipy import stats as st\n",
    "from tkinter import *\n",
    "from matplotlib.figure import Figure\n",
    "from matplotlib.backends.backend_tkagg import (FigureCanvasTkAgg, NavigationToolbar2Tk)\n",
    "import warnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9bb5bb8e-ba21-46a7-9470-1ab0e3cd67ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "from image import *\n",
    "from process import *\n",
    "from classes import *\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "048a1d09-912d-41e9-8edd-456e8895b5ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Debugging \n",
    "img_dict = {}\n",
    "width_dict = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0d1058ba-dd86-4c8d-8d8f-28adaef549e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "warnings.filterwarnings(\"ignore\")\n",
    "device = torch.device(\"cpu\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8c5b7586-5c90-41f8-8dba-8b656f121aef",
   "metadata": {},
   "outputs": [],
   "source": [
    "green_filename = 'h5_greenonly20230111_4.h5' #green dataset\n",
    "red_filename = 'h5_redonly20230111_4.h5' # red dataset, use red for\n",
    "\n",
    "file_path_model = \"D:/Models/ver15.pth\"\n",
    "file_path_data = \"D:/Datasets/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0c9ea88f-309b-418e-990f-2e487e1eb3f0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SegNet(\n",
       "  (conv1_1): KerasConv(\n",
       "    (conv): Conv2d(1, 64, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "    (activation): ReLU()\n",
       "  )\n",
       "  (conv1_2): KerasConv(\n",
       "    (conv): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "    (activation): ReLU()\n",
       "  )\n",
       "  (pool1): MaxPool2d(kernel_size=(2, 2), stride=(2, 2), padding=0, dilation=1, ceil_mode=False)\n",
       "  (layer_1): Sequential(\n",
       "    (0): KerasConv(\n",
       "      (conv): Conv2d(1, 64, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "    (1): KerasConv(\n",
       "      (conv): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "  )\n",
       "  (conv2_1): KerasConv(\n",
       "    (conv): Conv2d(64, 128, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "    (activation): ReLU()\n",
       "  )\n",
       "  (conv2_2): KerasConv(\n",
       "    (conv): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "    (activation): ReLU()\n",
       "  )\n",
       "  (pool2): MaxPool2d(kernel_size=(2, 2), stride=(2, 2), padding=0, dilation=1, ceil_mode=False)\n",
       "  (layer_2): Sequential(\n",
       "    (0): KerasConv(\n",
       "      (conv): Conv2d(64, 128, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "    (1): KerasConv(\n",
       "      (conv): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "  )\n",
       "  (conv3_1): KerasConv(\n",
       "    (conv): Conv2d(128, 256, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "    (activation): ReLU()\n",
       "  )\n",
       "  (conv3_2): KerasConv(\n",
       "    (conv): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "    (activation): ReLU()\n",
       "  )\n",
       "  (pool3): MaxPool2d(kernel_size=(2, 2), stride=(2, 2), padding=0, dilation=1, ceil_mode=False)\n",
       "  (layer_3): Sequential(\n",
       "    (0): KerasConv(\n",
       "      (conv): Conv2d(128, 256, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "    (1): KerasConv(\n",
       "      (conv): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "  )\n",
       "  (conv4_1): KerasConv(\n",
       "    (conv): Conv2d(256, 512, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "    (activation): ReLU()\n",
       "  )\n",
       "  (conv4_2): KerasConv(\n",
       "    (conv): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "    (activation): ReLU()\n",
       "  )\n",
       "  (drop4): Dropout(p=0.5, inplace=False)\n",
       "  (pool4): MaxPool2d(kernel_size=(2, 2), stride=(2, 2), padding=0, dilation=1, ceil_mode=False)\n",
       "  (layer_4): Sequential(\n",
       "    (0): KerasConv(\n",
       "      (conv): Conv2d(256, 512, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "    (1): KerasConv(\n",
       "      (conv): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "    (2): Dropout(p=0.5, inplace=False)\n",
       "  )\n",
       "  (conv5_1): KerasConv(\n",
       "    (conv): Conv2d(512, 1024, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "    (activation): ReLU()\n",
       "  )\n",
       "  (conv5_2): KerasConv(\n",
       "    (conv): Conv2d(1024, 1024, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "    (activation): ReLU()\n",
       "  )\n",
       "  (drop5): Dropout(p=0.5, inplace=False)\n",
       "  (layer_5): Sequential(\n",
       "    (0): KerasConv(\n",
       "      (conv): Conv2d(512, 1024, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "    (1): KerasConv(\n",
       "      (conv): Conv2d(1024, 1024, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "    (2): Dropout(p=0.5, inplace=False)\n",
       "  )\n",
       "  (layer_6): KerasUp(\n",
       "    (upsample): Upsample(scale_factor=(2.0, 2.0), mode=nearest)\n",
       "    (conv1): KerasConv(\n",
       "      (conv): Conv2d(1024, 512, kernel_size=(2, 2), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "    (activation): ReLU()\n",
       "    (conv2): KerasConv(\n",
       "      (conv): Conv2d(1024, 512, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "    (conv3): KerasConv(\n",
       "      (conv): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "  )\n",
       "  (layer_7): KerasUp(\n",
       "    (upsample): Upsample(scale_factor=(2.0, 2.0), mode=nearest)\n",
       "    (conv1): KerasConv(\n",
       "      (conv): Conv2d(512, 256, kernel_size=(2, 2), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "    (activation): ReLU()\n",
       "    (conv2): KerasConv(\n",
       "      (conv): Conv2d(512, 256, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "    (conv3): KerasConv(\n",
       "      (conv): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "  )\n",
       "  (layer_8): KerasUp(\n",
       "    (upsample): Upsample(scale_factor=(2.0, 2.0), mode=nearest)\n",
       "    (conv1): KerasConv(\n",
       "      (conv): Conv2d(256, 128, kernel_size=(2, 2), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "    (activation): ReLU()\n",
       "    (conv2): KerasConv(\n",
       "      (conv): Conv2d(256, 128, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "    (conv3): KerasConv(\n",
       "      (conv): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "  )\n",
       "  (layer_9): KerasUp(\n",
       "    (upsample): Upsample(scale_factor=(2.0, 2.0), mode=nearest)\n",
       "    (conv1): KerasConv(\n",
       "      (conv): Conv2d(128, 64, kernel_size=(2, 2), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "    (activation): ReLU()\n",
       "    (conv2): KerasConv(\n",
       "      (conv): Conv2d(128, 64, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "    (conv3): KerasConv(\n",
       "      (conv): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "  )\n",
       "  (layer_10): Sequential(\n",
       "    (0): KerasConv(\n",
       "      (conv): Conv2d(64, 2, kernel_size=(3, 3), stride=(1, 1), padding=same)\n",
       "      (activation): ReLU()\n",
       "    )\n",
       "    (1): Conv2d(2, 1, kernel_size=(1, 1), stride=(1, 1), padding=same)\n",
       "  )\n",
       ")"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Set up segmentation model\n",
    "segnet = SegNet()\n",
    "segnet.load_state_dict(torch.load(file_path_model,map_location=device))\n",
    "segnet.eval()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "6fc310e6-5678-4105-86e1-4dabf018d667",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Load Data\n",
    "r = h5py.File(file_path_data + red_filename, 'r') # Load h5\n",
    "g = h5py.File(file_path_data + green_filename,\"r\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f88cf914-bc1c-4aee-80c0-3b344b434de1",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_frames = len(r.keys())-1\n",
    "spacing = 3\n",
    "start = 5300\n",
    "j = '05300'\n",
    "showUI = True\n",
    "\n",
    "activation = {}\n",
    "\n",
    "if not showUI:\n",
    "    first_intestine_ind = 1310"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2b932774-9bf0-46c7-917a-dbc5d0936863",
   "metadata": {},
   "outputs": [],
   "source": [
    "lag_lengths = np.zeros(total_frames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "56582bdf-352e-42ea-99bd-4cd049d91766",
   "metadata": {},
   "outputs": [],
   "source": [
    "#UI functions \n",
    "\n",
    "canvas = None\n",
    "fig = Figure()\n",
    "window = Tk()\n",
    "window.title(\"Identify pharyngeal intestine junction\")\n",
    "\n",
    "\n",
    "\n",
    "def plotGraphUI():\n",
    "    global canvas, ROI\n",
    "    canvas = FigureCanvasTkAgg(fig, master = window)\n",
    "\n",
    "    plot1 = fig.add_subplot(111)\n",
    "    flatten = ROI.flatten()\n",
    "    #mode = st.mode(np.rint(flatten[flatten > 20]))[0][0] + threshold\n",
    "    #print(mode)\n",
    "    FI = ROI.copy()\n",
    "    #FI[FI < mode] = 0 #try no threshold\n",
    "    plot1.imshow(FI)\n",
    "    #To get xticks range\n",
    "    max_x = FI.shape[1]\n",
    "    plot1.set_xticks(range(0,max_x, 100))\n",
    "\n",
    "    canvas.draw()\n",
    "    canvas.get_tk_widget().pack(fill = BOTH, expand = True)\n",
    "\n",
    "def clearGraphUI():\n",
    "    global canvas\n",
    "    print(\"cleared\")\n",
    "    canvas.get_tk_widget().destroy()\n",
    "    cavas = None\n",
    "\n",
    "\n",
    "def getInput():\n",
    "    global plot_exists, first_intestine_ind\n",
    "    input = textbox.get(1.0, \"end-1c\")\n",
    "    first_intestine_ind = int(input)\n",
    "    #print(plot_exists)\n",
    "    if plot_exists:\n",
    "        clearGraphUI()\n",
    "    plotGraphUI()\n",
    "    plot_exists = True\n",
    "    \n",
    "textbox = Text(window, height = 5, width = 20)\n",
    "textbox.pack()\n",
    "plot_exists = False\n",
    "plot_btn = Button(window, command = getInput, text = \"Identify intestine region\")\n",
    "plot_btn.pack()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "18bf896e-71ae-470b-a589-41310288d59f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.2486902203969274\n"
     ]
    }
   ],
   "source": [
    "#Process first frame\n",
    "arr = get_image(r, 0)\n",
    "green_arr = get_image(g,0)\n",
    "mask = getPrediction(arr, segnet, torch.device(\"cpu\"))\n",
    "mask = ChooseLargestBlob(mask)\n",
    "mask_filled = fillHoles(mask)\n",
    "mask_erode = erode(mask_filled)\n",
    "mask_skeleton = skeletonize(mask_erode)\n",
    "body_points = sortSkeleton(mask_skeleton)\n",
    "original_img = padEdges(arr)\n",
    "original_green_img = padEdges(green_arr)\n",
    "\n",
    "[spline_x, spline_y, t_vals, x_vals, y_vals] = fitSpline(body_points)\n",
    "ROI = displayThreshold(original_img,  spline_x, spline_y, t_vals, x_vals, y_vals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "93dc33e1-1868-4ef6-a64a-eb3cffb436c3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cleared\n"
     ]
    }
   ],
   "source": [
    "if showUI:\n",
    "    window.mainloop()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "8ef1d04c-8055-40c5-a656-e70239cee0bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "[t_vals_1, body_points] = pointsOnSpline(spline_x, spline_y, max(t_vals), resolution = 4)\n",
    "arc_lengths = arcLengthSpline(spline_x, spline_y, spline_x(t_vals_1), spline_y(t_vals_1)) #t_vals is 4 times longer than x_vals\n",
    "[integrated_vals, ROI_l, width_dist] = integrateImage([original_img, original_green_img], spline_x, spline_y, t_vals_1, arc_lengths)\n",
    "xpoints_1 = arc_lengths #variable switch here because I was too lazy\n",
    "total_val = integrated_vals[:,0]\n",
    "total_green_val = integrated_vals[:,1]\n",
    "\n",
    "max_arc_length = math.floor(0.9*xpoints_1.max()) #first frame sets the kymograph y dim for all frames.\n",
    "[waveform,green_waveform] = interpolateIntensity([total_val, total_green_val], xpoints_1, max_arc_length, spacing)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "43cc770a-6338-4929-9152-5f21cd0385e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "img_dict[start] = ROI\n",
    "red_kym = np.zeros([total_frames, max_arc_length*spacing+400])\n",
    "green_kym = np.zeros([total_frames, max_arc_length*spacing+400])\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "8f227c61-bf6b-464b-a5d6-98553b695f4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "intestine_target = int(xpoints_1[first_intestine_ind])\n",
    "\n",
    "kymo_length = max_arc_length *spacing\n",
    "        #plotting just entry is good enough\n",
    "\n",
    "#For alignment\n",
    "kernel7 = np.tile(np.arange(-3,4),(7,1))\n",
    "res7_1 = signal.correlate2d(ROI_l, kernel7)\n",
    "avg_act1 = np.mean(res7_1, axis = 0)\n",
    "lag_lengths[start] = np.argmax(avg_act1[intestine_target-20:intestine_target+20]) + intestine_target - 20\n",
    "#plt.imshow(ROI)\n",
    "#plt.show()\n",
    "#plt.plot(avg_act1)\n",
    "#plt.show()\n",
    "red_kym[start] = waveform\n",
    "green_kym[start] = green_waveform\n",
    "width_dict[start] = width_dist\n",
    "activation[start] = avg_act1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "17cc7ef5-9cbc-49fa-b46b-9892c6b9dd42",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "global var?  (2836,) (2528,)\n",
      "(2836, 2) (2528,)\n",
      "(2528,) (2836,) (2836,)\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "x and y arrays must be equal in length along interpolation axis.",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "Input \u001b[1;32mIn [34]\u001b[0m, in \u001b[0;36m<cell line: 1>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     71\u001b[0m xpoints \u001b[38;5;241m=\u001b[39m xpoints \u001b[38;5;241m+\u001b[39m (intestine_target \u001b[38;5;241m-\u001b[39m intestine2_len)\n\u001b[0;32m     73\u001b[0m \u001b[38;5;28mprint\u001b[39m(xpoints\u001b[38;5;241m.\u001b[39mshape, total_val\u001b[38;5;241m.\u001b[39mshape, total_green_val\u001b[38;5;241m.\u001b[39mshape)\n\u001b[1;32m---> 74\u001b[0m [waveform,green_waveform] \u001b[38;5;241m=\u001b[39m \u001b[43minterpolateIntensity\u001b[49m\u001b[43m(\u001b[49m\u001b[43m[\u001b[49m\u001b[43mtotal_val\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtotal_green_val\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mxpoints\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmax_arc_length\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mspacing\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m     76\u001b[0m red_kym[k] \u001b[38;5;241m=\u001b[39m waveform\n\u001b[0;32m     77\u001b[0m green_kym[k] \u001b[38;5;241m=\u001b[39m green_waveform\n",
      "File \u001b[1;32m~\\Desktop\\Zhen Lab\\2022\\Kymograph_v1\\Kymograph_v1\\process.py:325\u001b[0m, in \u001b[0;36minterpolateIntensity\u001b[1;34m(values, xpoints, max_arc_length, spacing)\u001b[0m\n\u001b[0;32m    323\u001b[0m dict_func \u001b[38;5;241m=\u001b[39m {}\n\u001b[0;32m    324\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m val \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(\u001b[38;5;28mlen\u001b[39m(values)):\n\u001b[1;32m--> 325\u001b[0m     dict_func[val] \u001b[38;5;241m=\u001b[39m \u001b[43minterp1d\u001b[49m\u001b[43m(\u001b[49m\u001b[43mxpoints\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mvalues\u001b[49m\u001b[43m[\u001b[49m\u001b[43mval\u001b[49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    327\u001b[0m \u001b[38;5;124;03m'''Change based on expected length'''\u001b[39;00m\n\u001b[0;32m    328\u001b[0m xnew \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mlinspace(\u001b[38;5;241m1\u001b[39m,max_arc_length,max_arc_length\u001b[38;5;241m*\u001b[39mspacing)\n",
      "File \u001b[1;32m~\\Desktop\\virtualenvs\\research\\lib\\site-packages\\scipy\\interpolate\\_interpolate.py:446\u001b[0m, in \u001b[0;36minterp1d.__init__\u001b[1;34m(self, x, y, kind, axis, copy, bounds_error, fill_value, assume_sorted)\u001b[0m\n\u001b[0;32m    442\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__init__\u001b[39m(\u001b[38;5;28mself\u001b[39m, x, y, kind\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mlinear\u001b[39m\u001b[38;5;124m'\u001b[39m, axis\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m,\n\u001b[0;32m    443\u001b[0m              copy\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m, bounds_error\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m, fill_value\u001b[38;5;241m=\u001b[39mnp\u001b[38;5;241m.\u001b[39mnan,\n\u001b[0;32m    444\u001b[0m              assume_sorted\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m):\n\u001b[0;32m    445\u001b[0m     \u001b[38;5;124;03m\"\"\" Initialize a 1-D linear interpolation class.\"\"\"\u001b[39;00m\n\u001b[1;32m--> 446\u001b[0m     \u001b[43m_Interpolator1D\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[38;5;21;43m__init__\u001b[39;49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43my\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43maxis\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43maxis\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    448\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mbounds_error \u001b[38;5;241m=\u001b[39m bounds_error  \u001b[38;5;66;03m# used by fill_value setter\u001b[39;00m\n\u001b[0;32m    449\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcopy \u001b[38;5;241m=\u001b[39m copy\n",
      "File \u001b[1;32m~\\Desktop\\virtualenvs\\research\\lib\\site-packages\\scipy\\interpolate\\_polyint.py:54\u001b[0m, in \u001b[0;36m_Interpolator1D.__init__\u001b[1;34m(self, xi, yi, axis)\u001b[0m\n\u001b[0;32m     52\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdtype \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m\n\u001b[0;32m     53\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m yi \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m---> 54\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_set_yi\u001b[49m\u001b[43m(\u001b[49m\u001b[43myi\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mxi\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mxi\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43maxis\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43maxis\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32m~\\Desktop\\virtualenvs\\research\\lib\\site-packages\\scipy\\interpolate\\_polyint.py:124\u001b[0m, in \u001b[0;36m_Interpolator1D._set_yi\u001b[1;34m(self, yi, xi, axis)\u001b[0m\n\u001b[0;32m    122\u001b[0m     shape \u001b[38;5;241m=\u001b[39m (\u001b[38;5;241m1\u001b[39m,)\n\u001b[0;32m    123\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m xi \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m \u001b[38;5;129;01mand\u001b[39;00m shape[axis] \u001b[38;5;241m!=\u001b[39m \u001b[38;5;28mlen\u001b[39m(xi):\n\u001b[1;32m--> 124\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mx and y arrays must be equal in length along \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m    125\u001b[0m                      \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124minterpolation axis.\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m    127\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_y_axis \u001b[38;5;241m=\u001b[39m (axis \u001b[38;5;241m%\u001b[39m yi\u001b[38;5;241m.\u001b[39mndim)\n\u001b[0;32m    128\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_y_extra_shape \u001b[38;5;241m=\u001b[39m yi\u001b[38;5;241m.\u001b[39mshape[:\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_y_axis]\u001b[38;5;241m+\u001b[39myi\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_y_axis\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m:]\n",
      "\u001b[1;31mValueError\u001b[0m: x and y arrays must be equal in length along interpolation axis."
     ]
    }
   ],
   "source": [
    "\n",
    "for j in range(start+1, 5500): #Loop over how many frames you want to analyse\n",
    "    k = j\n",
    "    arr = np.zeros(np.shape(r['t00000']['s00']['0']['cells'])[1:])\n",
    "    green_arr = np.zeros(np.shape(g['t00000']['s00']['0']['cells'])[1:])\n",
    "    j = str(j)\n",
    "    lenj = len(j)\n",
    "    for instance in range(5-lenj):\n",
    "        j = '0' + j\n",
    "    #print(j)\n",
    "    for i in range(np.shape(r['t00000']['s00']['0']['cells'])[0]): #Sum up z stack\n",
    "        t = r['t'+ j]['s00']['0']['cells']\n",
    "        t_g = g['t'+ j]['s00']['0']['cells']\n",
    "        green_arr = green_arr + t_g[i]\n",
    "        arr = arr + t[i] # arr holds the summed up z stack image\n",
    "        #print(arr)\n",
    "\n",
    "    start_time = time.time()\n",
    "    arr = arr.T #for test only#start_time = time.time()\n",
    "    green_arr = green_arr.T\n",
    "    original_dim = arr.shape\n",
    "\n",
    "##Over here we select whether it is red or green we integrate\n",
    "    original_arr = np.array(arr) #red\n",
    "    original_green_arr = np.array(green_arr) #green\n",
    "\n",
    "    mask = getPrediction(original_arr, segnet, torch.device(\"cpu\"))\n",
    "    mask = ChooseLargestBlob(mask)\n",
    "    mask_filled = fillHoles(mask)\n",
    "    mask_erode = erode(mask_filled)\n",
    "    mask_skeleton = skeletonize(mask_erode)\n",
    "    body_points = sortSkeleton(mask_skeleton)\n",
    "    original_img = padEdges(original_arr)\n",
    "    original_green_img = padEdges(original_green_arr)\n",
    "\n",
    "    [spline_x, spline_y, t_vals, x_vals, y_vals] = fitSpline(body_points)\n",
    "\n",
    "    [t_vals, body_points] = pointsOnSpline(spline_x, spline_y, max(t_vals), resolution = 4)\n",
    "    x_points = arcLengthSpline(spline_x, spline_y, spline_x(t_vals), spline_y(t_vals)) #t_vals is 4 times longer than x_vals\n",
    "    #print(\"global var? \",arcLengthSpline(spline_x, spline_y, spline_x(t_vals), spline_y(t_vals)).shape, xpoints.shape)\n",
    "    [integrated_vals, ROI_l, width_dist] = integrateImage([original_img, original_green_img], spline_x, spline_y, t_vals, x_points)\n",
    "    #print(integrated_vals.shape, xpoints.shape)\n",
    "    cur_arc_length = body_points.max()\n",
    "    \n",
    "    total_val = integrated_vals[:,0]\n",
    "    total_green_val = integrated_vals[:,1]\n",
    "\n",
    "\n",
    "    #Cross correlate entire waveform to get approximate displacement vector\n",
    "    res7 = signal.correlate2d(ROI_l, kernel7)\n",
    "    avg_act = np.mean(res7, axis = 0)\n",
    "    corr = signal.correlate(avg_act[20:-20], activation[k-1][20:-20], mode = 'full')\n",
    "    lags = signal.correlation_lags(len(avg_act[20:-20]),len(activation[k-1][20:-20]), mode='full')\n",
    "    lag = lags[np.argmax(corr)]\n",
    "    #print(lag) #Should be 19\n",
    "    smooth_actROI = signal.savgol_filter(avg_act, 15, 3)[int(lag_lengths[k-1])+lag-40:int(lag_lengths[k-1])+lag+40]\n",
    "    maxima = np.array(signal.argrelextrema(smooth_actROI, np.greater)[0].tolist())\n",
    "\n",
    "    prev_loc = np.array([40,activation[k-1][int(lag_lengths[k-1])]/1000])\n",
    "\n",
    "    maximaXandY = (np.concatenate((np.expand_dims(maxima,axis = 1),np.expand_dims(smooth_actROI[maxima]/1000,axis = 1)),axis = 1))\n",
    "    maximaDist = np.linalg.norm(maximaXandY - prev_loc,axis = 1) #y is distorted, divide by 100\n",
    "    intestine2_len = maxima[np.argmin(maximaDist)] + int(lag_lengths[k-1])+lag-40 # Just x distance doesnt work, use y distance\n",
    "\n",
    "    activation[k] = avg_act\n",
    "\n",
    "\n",
    "\n",
    "    lag_lengths[k] = intestine2_len\n",
    "\n",
    "    #Need to factor in previous start point\n",
    "    xpoints = xpoints + (intestine_target - intestine2_len)\n",
    "    \n",
    "    print(xpoints.shape, total_val.shape, total_green_val.shape)\n",
    "    [waveform,green_waveform] = interpolateIntensity([total_val, total_green_val], xpoints, max_arc_length, spacing)\n",
    "\n",
    "    red_kym[k] = waveform\n",
    "    green_kym[k] = green_waveform\n",
    "    secs = time.gmtime((total_frames - k)*(time.time()-start_time))\n",
    "    time_remaining = time.strftime(\"%H:%M:%S\", secs)\n",
    "    print(\"Done t=\" + str(k) + \" out of \"+ str(total_frames) + \"Estimated time remaining:\" + str(time_remaining))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "559ddf37-8aaf-418e-80e1-74baf0e9527b",
   "metadata": {},
   "outputs": [],
   "source": [
    "red_kym = red_kym.T\n",
    "green_kym = green_kym.T\n",
    "r.close()\n",
    "g.close()\n",
    "\n",
    "np.save(\"red_kym.npy\", red_kym)\n",
    "np.save(\"green_kym.npy\", green_kym)\n",
    "np.save(\"lag_indices.npy\", lag_lengths)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "research",
   "language": "python",
   "name": "research"
  },
  "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.10.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}