4_refine_FOV_to_exvivo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4a73fcac",
   "metadata": {},
   "source": [
    "# refine FOV to exvivo registration and save the matching cells coordinates\n",
    "- requires napari to refine registration and selecting the dots to save."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "65f42e01",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "070924\n"
     ]
    }
   ],
   "source": [
    "\n",
    "from datetime import date\n",
    "import tifffile, pickle, napari, pprint\n",
    "import sys, glob\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "np.set_printoptions(suppress=True)\n",
    "import scipy.ndimage as ndi\n",
    "\n",
    "from datetime import date\n",
    "today = date.today()\n",
    "todaystamp = today.strftime(\"%m%d%y\")\n",
    "print(todaystamp)\n",
    "sys.path.append('./functions/')\n",
    "import transformations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0750ca4",
   "metadata": {},
   "source": [
    "# load exvivo "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c7159701",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 48, 1517, 1618)\n"
     ]
    }
   ],
   "source": [
    "datadir= '/path_to_the_data/'\n",
    "exvivo_f = datadir + 'exvivo_2xds.pkl'\n",
    "exvivo = pickle.load(open(exvivo_f,'rb'))[:,::-1]  #144um, because each z is 1.5 um!\n",
    "exvivo[exvivo<=0]=0\n",
    "print(exvivo.shape)  # channels are BV/tomatos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a634d390",
   "metadata": {},
   "outputs": [],
   "source": [
    "Z = exvivo.shape[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "008a8ea8",
   "metadata": {},
   "source": [
    "# also load transformed invivo as the reference \n",
    "- note here we are loading the output from the transformed result from `invivo_to_exvivo` notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a39d7e28",
   "metadata": {},
   "outputs": [],
   "source": [
    "base_dir='/set_your_output_directory/'  # note this is also where the output will be stored!\n",
    "invivo_t = pickle.load(open(base_dir + f'transformed_06202024.pkl','rb'))['transformed'][:Z]\n",
    "# invivo_t = pickle.load(open('/Users/shuonan/Downloads/invivo_tomato_T_2xds.pkl', 'rb'))[:Z]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "af1607d2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(48, 1517, 1618)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "invivo_t.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c077364",
   "metadata": {},
   "source": [
    "# load ONE fov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f29ebb32",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['/Users/shuonan/Downloads/transformed_FOVs/FOV_Stacked_ch2_movie_FOV13_14_zstack_1100nm_3x_.tif_mapped_to_exvivo_2xds.pkl',\n",
       " '/Users/shuonan/Downloads/transformed_FOVs/FOV_Stacked_ch2_movie_FOV9_zstack_1100nm_3x_.tif_mapped_to_exvivo_2xds.pkl',\n",
       " '/Users/shuonan/Downloads/transformed_FOVs/FOV_Stacked_ch2_movie_FOV3_4_zstack_1100nm_3x_.tif_mapped_to_exvivo_2xds.pkl',\n",
       " '/Users/shuonan/Downloads/transformed_FOVs/FOV_Stacked_ch2_movie_FOV11_12_zstack_1100nm_3x_.tif_mapped_to_exvivo_2xds.pkl',\n",
       " '/Users/shuonan/Downloads/transformed_FOVs/FOV_Stacked_ch2_movie_FOV10_zstack_1100nm_3x.tif_mapped_to_exvivo_2xds.pkl',\n",
       " '/Users/shuonan/Downloads/transformed_FOVs/FOV_Stacked_ch2_movie_FOV1_2_zstack_1100nm_3x_.tif_mapped_to_exvivo_2xds.pkl']"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "transformedFOV_path = \"/Users/shuonan/Downloads/transformed_FOVs/\"\n",
    "glob.glob(transformedFOV_path+\"*\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "de30c01e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(48, 1517, 1618)\n"
     ]
    }
   ],
   "source": [
    "fovfile = transformedFOV_path + 'FOV_Stacked_ch2_movie_FOV11_12_zstack_1100nm_3x_.tif_mapped_to_exvivo_2xds.pkl'\n",
    "mapped_fov = pickle.load(open(fovfile, 'rb'))[:Z]\n",
    "print(mapped_fov.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "110c6c85",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fdc64817a60>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAisAAAA/CAYAAAAhdoiMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAU60lEQVR4nO3df1AU5/0H8Pezu/cD8LgI9ICDiPiNViPGGLTaxkaTOsbGVI2d1lhDdDrTqR01knSsWttJ2mmD/adN24lkkmacydh8ZTrRxCZOEmyV6EDCd0AqamMwQUEjXiJw/L67vf18/9i75VbwFzlwOT6vmRth97m9530My8dnn2dPEBGBMcYYY8yipNvdAcYYY4yx6+FihTHGGGOWxsUKY4wxxiyNixXGGGOMWRoXK4wxxhizNC5WGGOMMWZpXKwwxhhjzNK4WGGMMcaYpXGxwhhjjDFL42KFMcYYY5Y2bMXKrl27kJ+fD6fTicLCQhw9enS4XooxxhhjCWxYipWysjIUFxdjx44dOH78OL797W/ju9/9Lpqamobj5RhjjDGWwMRwfJDh3Llzcd9996G0tNTYNm3aNKxYsQIlJSXxfjnGGGOMJTAl3gcMBoOoqanBtm3bTNsXL16MysrKAe0DgQACgYDxvaZpaG1tRXp6OoQQ8e4eY4wxxoYBEaGzsxNerxeSFN8LN3EvVr788kuEw2FkZmaatmdmZqKlpWVA+5KSEvzmN7+JdzcYY4wxdhs0NzcjNzc3rseMe7ESdfWoCBENOlKyfft2PPPMM8b3fr8fEyZMwHw8AgW24eoeYwlFycuF5kqGCOtXdSn2d00WgBD6NhkgWdIfkt6mfYoTqY1BOHzd6JzsxrjGTojLrUA4rD9fSEBqCiBJQEcXhCSBNE3/fZYlaOlu9HpTkNTcBSgSpI5uvQ+y3N+HSH+0FCeCaQ44LvdA9AWANr+xD1oYwuEEBfoASUbborvguz+MKX/rhrjcCiEEwldaQao6vG8mY2xIVIRwDAfhcrnifuy4FysZGRmQZXnAKIrP5xsw2gIADocDDodjkI7ZoAguVhi7GYqwg2AHpEixYpMhVA0kC0CSEE6yAYoEUgTUFAWCCEIlyAENKd02yE4Ztq52pNW36c+X7ACFgeiUto4+vaiQ7RA2W38hoyhAlwrbZ11om5OJ1ukC//O/rUBI7S9CZKm/WJEVyH0yJLsTAgqgBICwBkhCfz1hAxS9ref/WpFxJhniiy4ANkAICGEzF2KMMeuInC6GYwpH3IsVu92OwsJClJeX47HHHjO2l5eXY/ny5fF+OcYY0F8YhDUAACXb9T/uQRUIq1ACIUCSoKY60TPRAVezPpIiVA2Oc316URJ5aO4UiO7e/uOGNYA0fYQFABQFFArp30dGWIQQcLaFITQFCIYggqH+PkWuXZMkIAWCQLf+vQipIKJI4RMzCqMRQBqoswuiqzummBEgLe7rARhjo8AtzYApKSnBnDlz4HK54PF4sGLFCpw5c8bUZt26daiurkZpaalxEsvOzkZTUxPWr18f184zxiJCqv4Ih4FwGCQJhMY7ATUMEd0XUqF09MF1MQQ1RYaa6tQLgUBQf4TDoL4ApNZOvTiJFjCAXkAA+vdC9BcuUURIPvk5Ju3xQfQF9fbhmGNoGkRYg9Aiz9e0yCN63MjXqp6BwpqeJXrJR4vpC2NszLmlYqWiogIbNmzAhx9+iPLycqiqisWLF6O7u9vUbsmSJXj++eeRm5sLm80Gj8eDgwcPIi8vL66dZ4zFiBl6lbsCUDqD+lwVOebXXNPgvNiJULKElnnJ+iUa0yEEKBCMtCX9IQlAlgFJQNjt6CzMAWWmG+31y0P6yIjoDUTaS4AiA7Ksz49RZOMBWdLby5GvpejlIhmQ5MHnt0nClI8xNrbc0mWgd9991/T97t274fF4UFNTgwceeMDY7nA4sH37dmzfvv2Gx7x66XJHR8etdIkxFiVL/SMfRJA6+/Tt0T/yFBntUMNwn+2Gs9UBBPXLQ9A0Y+RCyBKIRGREhPpHPwDAYUfbZAVS0IXkli9AYQB3uPTRm0BQPxYJ0+uKmIKIFBlQ0f96QgIcDqC3D0DYmAtDRBBa9FmRy0T84SCMjVlf6dff7/cDANLS0kzbjxw5Ao/HgylTpuAnP/kJfD7fNY9RUlICt9ttPO68886v0iXGxiTq64Po7oXoDUD0Bfr/jX4dDEFE5rNAkiD1BJHU5AckCRQMgYJBkKqCQiF9Pkow1D+JliJzVsJhUFc3co50IPnslUgxEwalOEFJDpAa1lfqDHiE9UdIhdBI70t029fGIzQ1Rx+BIQKFw8YlIlJVUDjcP68l2h/G2Jgz5DvYEhGWL1+OtrY20+f+lJWVYdy4ccjLy0NjYyN+/etfQ1VV1NTUDLrqZ7CRlTvvvBMLsZxXAzF2k+Q73BBJSdfYqV9uIZsC2BT90pAiARogwmHAd0Wf0yJF5qLEFgWS6B9ZiblcE30tCgQhHHZ9ezhsTKY1RnOEMF3C0VKTAQ2QevRJvcHcNPR67Eg90gDq6z8P9L9+ZDJu5NKR5u/gpcuMWZRKIRzBW/D7/UhNTY3rsYe8Gmjjxo04ceIEjh07Ztq+atUq4+uCggLMnj0beXl5eOedd7By5coBx7nW0mXG2C0QEuB0DJyESmQuIETM3A8J+hUWJXIakOT+UZRYknlCrbDboeamQ1Mk2D9t6S9SonNQYgsVIYz7uUAIdHzdDc0mML5OX95s/+wy7Bft+oiK8XqRAsXoZ+T1JXngxF7G2JgwpGJl06ZNOHDgAD744IMb3qUuOzsbeXl5aGhoGFIHGWM3IXKjNxG7agcw5q8AgAhrIEkvRARF/uirkSJB0ifQIgy9IIgtWDQC7JJ+aQgAFAV9X3NCDmggNbJiRxL6MQBjBMaYfDs+FaKjWx9JcUlo/zpg67kD405fAfX0AoGAebVRZD6LcT8VkiAEDV5IMcbGhFsqVogImzZtwv79+3HkyBHk5+ff8DlXrlxBc3MzsrOzh9xJxtgNEEEEQvofe1nSi5BosRCddAvohcDVn9kRjt78TdYn4EbniSBScEgCQlEic0g0oLcXKZ/5ITq6zZdkYu/FEt00Lhkd08Zj3Hk75EutSDvdDXtnEpQeTb/PiixBjEsB+Tv1fsT2j/qXOZMkxUy4ZYyNNbdUrGzYsAGvv/463nrrLbhcLuMutW63G0lJSejq6sJzzz2H73//+8jOzsa5c+fwy1/+EhkZGaYbxDHGhkH0nibBkHlkBYis0iGIyNcU1gsZoUbuaaIoEClJoHa/UajEFgkUCkUOJ/Tjf+7Tb1YpxV7ykQaOqhDB7leNy1EkBFIbOiG6+4C+yN1rI6+PYHBgIRXtezQi3xSOsTHploqV0tJSAMDChQtN23fv3o1169ZBlmXU19fjtddeQ3t7O7Kzs/Hggw+irKzspj8rIHqiVBEybt3LGLs+ohCEZIvc7RWmuSvksOt3lCUAJIMgAxoASYIgVb80pAUhelQQVEDSIpd0hLEKCGrIPNFWRC/bQC9ShIAQmv4c6Mem5CSI3i7YznSDoK9YxhetCIxPRu/XXXCd6gN6g0AfQdgUaBQcWGQJCSKmWNGECiKeYMuYFanQ/1MzxHU71zXk1UDD5cKFC7x8mTHGGBulhuNTly1XrGiahjNnzuDuu+9Gc3Nz3Jc/WVV0yTZnTnxjMTdn5syJbCzmHiwzEaGzsxNerxfSYJd0v4K4f5DhVyVJEnJycgAAqampY+YHH8WZx46xmJszjw1jMTMwNnNfndntdg/L6/BNCxhjjDFmaVysMMYYY8zSLFmsOBwOPPvss2PqzraceewYi7k589gwFjMDYzP3SGe23ARbxhhjjLFYlhxZYYwxxhiL4mKFMcYYY5bGxQpjjDHGLI2LFcYYY4xZGhcrjDHGGLM0yxUru3btQn5+PpxOJwoLC3H06NHb3aUhKykpwZw5c+ByueDxeLBixQqcOXPG1IaI8Nxzz8Hr9SIpKQkLFy7EqVOnTG0CgQA2bdqEjIwMpKSkYNmyZbhw4cJIRhmykpISCCFQXFxsbEvEzBcvXsQTTzyB9PR0JCcn495770VNTY2xPxEzq6qKX/3qV8jPz0dSUhImTZqE3/72t9A0zWgz2nN/8MEH+N73vgev1wshBN58803T/njla2trQ1FREdxuN9xuN4qKitDe3j7M6QZ3vcyhUAhbt27FjBkzkJKSAq/XiyeffBKff/656RijLTNw4591rJ/+9KcQQuCFF14wbR9tuW8m83//+18sW7YMbrcbLpcL8+bNQ1NTk7F/xDKThezdu5dsNhu98sordPr0adq8eTOlpKTQ+fPnb3fXhuThhx+m3bt308mTJ6muro6WLl1KEyZMoK6uLqPNzp07yeVy0RtvvEH19fW0atUqys7Opo6ODqPN+vXrKScnh8rLy6m2tpYefPBBmjlzJqmqejti3bTq6mqaOHEi3XPPPbR582Zje6Jlbm1tpby8PFq3bh199NFH1NjYSIcOHaKzZ88abRItMxHR7373O0pPT6e3336bGhsb6R//+AeNGzeOXnjhBaPNaM998OBB2rFjB73xxhsEgPbv32/aH698S5YsoYKCAqqsrKTKykoqKCigRx99dKRimlwvc3t7Oy1atIjKysro448/pqqqKpo7dy4VFhaajjHaMhPd+GcdtX//fpo5cyZ5vV7605/+ZNo32nLfKPPZs2cpLS2NtmzZQrW1tfTpp5/S22+/TZcvXzbajFRmSxUr3/jGN2j9+vWmbVOnTqVt27bdph7Fl8/nIwBUUVFBRESaplFWVhbt3LnTaNPX10dut5teeuklItJPDjabjfbu3Wu0uXjxIkmSRO++++7IBrgFnZ2dNHnyZCovL6cFCxYYxUoiZt66dSvNnz//mvsTMTMR0dKlS+nHP/6xadvKlSvpiSeeIKLEy331yTxe+U6fPk0A6MMPPzTaVFVVEQD6+OOPhznV9V3vj3ZUdXU1ATD+UznaMxNdO/eFCxcoJyeHTp48SXl5eaZiZbTnHizzqlWrjN/nwYxkZstcBgoGg6ipqcHixYtN2xcvXozKysrb1Kv48vv9AIC0tDQAQGNjI1paWkyZHQ4HFixYYGSuqalBKBQytfF6vSgoKLD0+7JhwwYsXboUixYtMm1PxMwHDhzA7Nmz8YMf/AAejwezZs3CK6+8YuxPxMwAMH/+fPzrX//CJ598AgD4z3/+g2PHjuGRRx4BkLi5o+KVr6qqCm63G3PnzjXazJs3D2632/LvAaCf14QQuOOOOwAkbmZN01BUVIQtW7Zg+vTpA/YnWm5N0/DOO+9gypQpePjhh+HxeDB37lzTpaKRzGyZYuXLL79EOBxGZmamaXtmZiZaWlpuU6/ih4jwzDPPYP78+SgoKAAAI9f1Mre0tMBut2P8+PHXbGM1e/fuRW1tLUpKSgbsS8TMn332GUpLSzF58mS89957WL9+PZ566im89tprABIzMwBs3boVq1evxtSpU2Gz2TBr1iwUFxdj9erVABI3d1S88rW0tMDj8Qw4vsfjsfx70NfXh23btuFHP/qR8cm7iZr5D3/4AxRFwVNPPTXo/kTL7fP50NXVhZ07d2LJkiV4//338dhjj2HlypWoqKgAMLKZla+QZVgIIUzfE9GAbaPRxo0bceLECRw7dmzAvqFktur70tzcjM2bN+P999+H0+m8ZrtEyqxpGmbPno3nn38eADBr1iycOnUKpaWlePLJJ412iZQZAMrKyrBnzx68/vrrmD59Ourq6lBcXAyv14u1a9ca7RIt99XikW+w9lZ/D0KhEB5//HFomoZdu3bdsP1ozlxTU4M///nPqK2tveX+jdbc0Ynyy5cvx9NPPw0AuPfee1FZWYmXXnoJCxYsuOZzhyOzZUZWMjIyIMvygErL5/MN+J/LaLNp0yYcOHAAhw8fRm5urrE9KysLAK6bOSsrC8FgEG1tbddsYyU1NTXw+XwoLCyEoihQFAUVFRX4y1/+AkVRjD4nUubs7Gzcfffdpm3Tpk0zZswn4s8ZALZs2YJt27bh8ccfx4wZM1BUVISnn37aGFFL1NxR8cqXlZWFy5cvDzj+F198Ydn3IBQK4Yc//CEaGxtRXl5ujKoAiZn56NGj8Pl8mDBhgnFeO3/+PH7+859j4sSJABIvd0ZGBhRFueG5baQyW6ZYsdvtKCwsRHl5uWl7eXk5vvWtb92mXn01RISNGzdi3759+Pe//438/HzT/vz8fGRlZZkyB4NBVFRUGJkLCwths9lMbS5duoSTJ09a8n35zne+g/r6etTV1RmP2bNnY82aNairq8OkSZMSLvP9998/YEn6J598gry8PACJ+XMGgJ6eHkiS+RQiy7LxP7JEzR0Vr3zf/OY34ff7UV1dbbT56KOP4Pf7LfkeRAuVhoYGHDp0COnp6ab9iZi5qKgIJ06cMJ3XvF4vtmzZgvfeew9A4uW22+2YM2fOdc9tI5r5pqfijoDo0uVXX32VTp8+TcXFxZSSkkLnzp273V0bkp/97GfkdrvpyJEjdOnSJePR09NjtNm5cye53W7at28f1dfX0+rVqwdd+pibm0uHDh2i2tpaeuihhyyztPNmxK4GIkq8zNXV1aQoCv3+97+nhoYG+vvf/07Jycm0Z88eo02iZSYiWrt2LeXk5BhLl/ft20cZGRn0i1/8wmgz2nN3dnbS8ePH6fjx4wSA/vjHP9Lx48eNlS/xyrdkyRK65557qKqqiqqqqmjGjBm3bTnr9TKHQiFatmwZ5ebmUl1dnem8FggEjGOMtsxEN/5ZX+3q1UBEoy/3jTLv27ePbDYbvfzyy9TQ0EB//etfSZZlOnr0qHGMkcpsqWKFiOjFF1+kvLw8stvtdN999xnLfEcjAIM+du/ebbTRNI2effZZysrKIofDQQ888ADV19ebjtPb20sbN26ktLQ0SkpKokcffZSamppGOM3QXV2sJGLmf/7zn1RQUEAOh4OmTp1KL7/8sml/Imbu6OigzZs304QJE8jpdNKkSZNox44dpj9aoz334cOHB/0dXrt2LRHFL9+VK1dozZo15HK5yOVy0Zo1a6itrW2EUppdL3NjY+M1z2uHDx82jjHaMhPd+Gd9tcGKldGW+2Yyv/rqq3TXXXeR0+mkmTNn0ptvvmk6xkhlFkRENz8OwxhjjDE2siwzZ4UxxhhjbDBcrDDGGGPM0rhYYYwxxpilcbHCGGOMMUvjYoUxxhhjlsbFCmOMMcYsjYsVxhhjjFkaFyuMMcYYszQuVhhjjDFmaVysMMYYY8zSuFhhjDHGmKX9PwHVQNmlTVZbAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(np.max(mapped_fov,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "f139e552",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11_12\n"
     ]
    }
   ],
   "source": [
    "fov_id = fovfile.split('/')[-1].split('FOV')[2].split('_zstack')[0]\n",
    "print(fov_id)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "837d349d",
   "metadata": {},
   "source": [
    "# launch napari \n",
    "- note we are NOT transforming invivo layer to save time. `in_layer` is purely for the diagnosis purpose, and only be used before starting the registration. feel free to remove it if necessary. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "989f2cca",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "learning scale..\n",
      "scale difference: 0.9815037085163179\n",
      "applying rigid..\n",
      "done transform!\n"
     ]
    }
   ],
   "source": [
    "## note here the transforming ones are the fov_layers!!!\n",
    "features = {'label': np.empty(0, dtype=int)}\n",
    "all_transform_FOV = []\n",
    "\n",
    "viewer = napari.view_image(exvivo[1][None],channel_axis=0,name=['EX Tomato'],\n",
    "                           contrast_limits=(0,exvivo[1].max()),colormap=\"green\", visible = [True], blending='additive')\n",
    "fov_layer = viewer.add_image(mapped_fov[None],channel_axis=0,name=[f'FOV {fov_id} tomato'],\n",
    "                             contrast_limits=(0,mapped_fov.max()), colormap='red',visible = [True], blending='additive')\n",
    "in_layer = viewer.add_image(invivo_t[None],channel_axis=0,name=[f'invivo tomato'],\n",
    "                             contrast_limits=(0,invivo_t.max()), colormap='red',visible = [True], blending='additive')\n",
    "\n",
    "\n",
    "pl_in = viewer.add_points(\n",
    "    size=10, edge_width=1, edge_color='red',face_color='transparent',name=f\"IN points\", text='label', \n",
    "    features=features,ndim=3,out_of_slice_display=True)\n",
    "@pl_in.events.data.connect\n",
    "def update_feature_default_invivo():  \n",
    "    global points_layer\n",
    "    no_of_points = len(pl_in.data)\n",
    "    pl_in.feature_defaults['label'] = no_of_points + 1\n",
    "    pl_in.properties[\"label\"][0:no_of_points] = range(1, no_of_points+1)\n",
    "    pl_in.text.values[0:no_of_points] = [str(i) for i in range(1, no_of_points+1)]\n",
    "    pl_in.text.color = 'red'\n",
    "    pl_in.text.translation = np.array([-10, 0])\n",
    "update_feature_default_invivo()\n",
    "pl_in.mode = 'add'\n",
    "\n",
    "pl_ex = viewer.add_points(\n",
    "    size=10, edge_width=1, edge_color='green',face_color='transparent',name=f\"EX points\", text='label', \n",
    "    features=features,ndim=3,out_of_slice_display=True)\n",
    "@pl_ex.events.data.connect\n",
    "def update_feature_default_exvivo():  \n",
    "    global points_layer\n",
    "    no_of_points = len(pl_ex.data)\n",
    "    pl_ex.feature_defaults['label'] = no_of_points + 1\n",
    "    pl_ex.properties[\"label\"][0:no_of_points] = range(1, no_of_points+1)\n",
    "    pl_ex.text.values[0:no_of_points] = [str(i) for i in range(1, no_of_points+1)]\n",
    "    pl_ex.text.color = 'green'\n",
    "    pl_ex.text.translation = np.array([-10, 0])\n",
    "update_feature_default_exvivo()\n",
    "pl_ex.mode = 'add'\n",
    "\n",
    "\n",
    "@viewer.bind_key('s', overwrite = True)\n",
    "def scale(viewer):  \n",
    "    C = 1\n",
    "    print('learning scale..')\n",
    "    m_invivo = viewer.layers[f\"IN points\"].data\n",
    "    m_exvivo = viewer.layers[f\"EX points\"].data\n",
    "    assert len(m_invivo)==len(m_exvivo), f\"{len(m_invivo), len(m_exvivo)}\"\n",
    "    s_invivo = np.sqrt(np.sum((m_invivo-np.mean(m_invivo,0))**2)/len(m_invivo))\n",
    "    s_exvivo = np.sqrt(np.sum((m_exvivo-np.mean(m_exvivo,0))**2)/len(m_exvivo))\n",
    "    scl = s_exvivo/s_invivo\n",
    "    print(f'scale difference: {scl}')\n",
    "    m_invivo_new = scl*m_invivo    \n",
    "    fov_layer_scaled = np.array([ndi.zoom(fov_layer[c].data, scl, order=1) for c in range(C)])    \n",
    "    for c in range(C):\n",
    "        fov_layer[c].data = fov_layer_scaled[c]\n",
    "    viewer.layers[\"IN points\"].data = m_invivo_new    \n",
    "    \n",
    "    \n",
    "@viewer.bind_key('t', overwrite = True)\n",
    "def transform(viewer):    \n",
    "    C=1\n",
    "    viewer.update_console(locals())\n",
    "    print('applying rigid..')\n",
    "    m_invivo = viewer.layers[f\"IN points\"].data\n",
    "    m_exvivo = viewer.layers[f\"EX points\"].data\n",
    "    assert len(m_invivo)==len(m_exvivo), f\"{len(m_invivo), len(m_exvivo)}\"\n",
    "    bhat = wahba(m_invivo, m_exvivo)\n",
    "    offset = -(bhat[:3,:3])@bhat[-1]\n",
    "    fov_layer_affined = np.array([ndi.affine_transform(fov_layer[c].data, bhat[:3,:3],\n",
    "                                                    output_shape = viewer.layers['EX Tomato'].data.shape,\n",
    "                                                    offset = offset, order=1) for c in range(C)])\n",
    "    foo = np.c_[m_invivo, np.ones((m_invivo.shape[0],1))]   #Nx4    \n",
    "    m_invivo_new = foo@bhat\n",
    "    viewer.layers[\"IN points\"].data = m_invivo_new    \n",
    "    all_transform_FOV.append(dict(bhat=bhat))\n",
    "    for c in range(C):\n",
    "        fov_layer[c].data = fov_layer_affined[c]\n",
    "    print(\"done transform!\")\n",
    "\n",
    "# again this is to deform in vivo points to map to exvivo points!\n",
    "@viewer.bind_key('d', overwrite = True)\n",
    "def deform(viewer):    \n",
    "    print('applying deformable..')\n",
    "    ksz = 10\n",
    "    m_invivo = viewer.layers[f\"IN points\"].data\n",
    "    m_exvivo = viewer.layers[f\"EX points\"].data \n",
    "    assert len(m_invivo)==len(m_exvivo), f\"{len(m_invivo), len(m_exvivo)}\"\n",
    "    shift = m_exvivo-m_invivo   # so the newcoords = old+shift --> vecfield will have how much shift needed to apply to the old coords\n",
    "    vec_field = np.zeros(fov_layer[0].data.shape + (3,))  # M1,M2,2\n",
    "    for p,loc in enumerate(m_invivo):\n",
    "        vec_field[int(loc[0]),int(loc[1]),int(loc[2])] = shift[p]\n",
    "    for c in range(3):\n",
    "        vec_field[...,c] = ndi.gaussian_filter(vec_field[...,c], ksz)\n",
    "    A = np.zeros_like(m_invivo)\n",
    "    for p,loc in enumerate(m_invivo):\n",
    "        A[p] = vec_field[int(loc[0]),int(loc[1]),int(loc[2])]\n",
    "    diag_step,_,_,_ = np.linalg.lstsq(A,shift,rcond=None)\n",
    "    step = np.diag(diag_step)\n",
    "    vec_field_total =vec_field*step;  # element wise. \n",
    "    all_transform_FOV.append(dict(vec_field_total=vec_field_total))\n",
    "    mapz_base, mapx_base, mapy_base = np.meshgrid(np.arange(fov_layer[0].data.shape[0]),np.arange(fov_layer[0].data.shape[1]), \n",
    "                                       np.arange(fov_layer[0].data.shape[2]), indexing='ij')\n",
    "    mapz=mapz_base-vec_field_total[:,:,:,0]\n",
    "    mapx=mapx_base-vec_field_total[:,:,:,1]\n",
    "    mapy=mapy_base-vec_field_total[:,:,:,2]\n",
    "    \n",
    "    for c in range(2):\n",
    "        print(f'... warping for channel {c}')\n",
    "        img_de = warp(fov_layer[c].data,np.array((mapz,mapx,mapy)), order = 1)            \n",
    "        fov_layer[c].data = img_de    \n",
    "    \n",
    "    m_invivo_new = np.zeros_like(m_invivo)  # POINTS\n",
    "    for p,loc in enumerate(m_invivo):\n",
    "        new_s = vec_field_total[int(loc[0]),int(loc[1]),int(loc[2])]\n",
    "        m_invivo_new[p] = loc+new_s\n",
    "    viewer.layers[\"IN points\"].data = m_invivo_new\n",
    "    \n",
    "    plt.figure(figsize = (10,3))\n",
    "    for c in range(3):\n",
    "        plt.subplot(1,3,c+1)\n",
    "        plt.imshow(np.max(vec_field_total[...,c], 0))\n",
    "        plt.colorbar()\n",
    "    print('done deformation!')\n",
    "    plt.show()\n",
    "        \n",
    "    \n",
    "@viewer.bind_key('c',overwrite = True)\n",
    "def clear_selected(viewer):        \n",
    "    viewer.layers[f\"IN points\"].data = np.empty((0, 3))\n",
    "    viewer.layers[f\"EX points\"].data = np.empty((0, 3))\n",
    "\n",
    "\n",
    "@viewer.bind_key('q', overwrite = True)\n",
    "def save_rez(viewer):     \n",
    "    out_dict = dict()\n",
    "    out_dict['transformed tomato']=viewer.layers[f'FOV {fov_id} tomato'].data.astype('float16')\n",
    "    out_dict['invivo points']=viewer.layers['IN points'].data\n",
    "    out_dict['exvivo points']=viewer.layers['EX points'].data\n",
    "#     out_dict['transformed gcamp']=viewer.layers[f'FOV {fov_id} gcamp'].data.astype('float16')\n",
    "    out_dict['transformations']=all_transform_FOV      \n",
    "    \n",
    "    pickle.dump(out_dict, open(base_dir + f'transformed_fov_{fov_id}_{todaystamp}.pkl','wb'))\n",
    "    print(f\"saved the results! ({base_dir + f'transformed_fov_{fov_id}_{todaystamp}.pkl'})\")\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73edbedb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7c01a21",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65947f4e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8010a17",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "cv2",
   "language": "python",
   "name": "cv2"
  },
  "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}