{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as patches\n",
    "from pathlib import Path\n",
    "import numpy as np\n",
    "import os\n",
    "import sys\n",
    "import glob\n",
    "import pandas as pd\n",
    "import xml.etree.ElementTree as et\n",
    "import datetime\n",
    "from skimage.io import imread, imsave\n",
    "from imageio import volread as imread\n",
    "\n",
    "import tifffile\n",
    "import pystackreg\n",
    "from pystackreg import StackReg\n",
    "from skimage.filters import threshold_otsu\n",
    "\n",
    "from pystackreg.util import to_uint16   # make sure version 0.2.5 (not anything below)\n",
    "#from ims_to_tiff import convert_to_tif  \n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "import seaborn as sns\n",
    "import pylab as pl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "CYCLE_NUMS = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Converting to TIF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cycles = glob.glob('Cycle_*') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cycles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "error_files = []\n",
    "for ims in cycles:\n",
    "    try:\n",
    "        sourceFile = ims\n",
    "        destFile = \"tif/\" + ims[:-4]+'.tif'\n",
    "        convert_to_tif(sourceFile, destFile)\n",
    "    except:\n",
    "        print(f'{ims} has B-tree error or import error')\n",
    "        error_files.append(ims)\n",
    "        pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "error_files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mkdir tif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NUM_FOVS = 225 - len(error_files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NUM_FOVS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove all erroneous files \n",
    "for i in error_files:\n",
    "    command = f'rm tif/*{i[-7:-4]}*'\n",
    "    ! {command}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for c in range(CYCLE_NUMS):\n",
    "    os.makedirs(f'tif/Cycle_{c}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for c in range(CYCLE_NUMS):\n",
    "    if c == 0:\n",
    "        command = 'mv tif/Cycle_F* tif/Cycle_0'\n",
    "    else:\n",
    "        command = f'mv tif/Cycle_{c}_F* tif/Cycle_{c}'    \n",
    "    ! {command}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Registration (Testing)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test image sizes\n",
    "im1 = imread('tif/Cycle_0/Cycle_F001.tif')\n",
    "im1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ref = imread('tif/Cycle_0/Cycle_F121.tif')\n",
    "mov = imread('tif/Cycle_1/Cycle_1_F121.tif')\n",
    " \n",
    "ref_max = ref.max(0)   # max proj of z for each channel in reference\n",
    "ref_binary = ref_max[0] > threshold_otsu(ref_max[0])  # binarizing for channels - T/F\n",
    "\n",
    "mov_max = mov.max(0)  # max projection by each cycle in for loop (ref above)\n",
    "mov_binary = mov_max[0] > threshold_otsu(mov_max[0]) # binary of moved image\n",
    "    \n",
    "sr = StackReg(StackReg.RIGID_BODY)\n",
    "tmat = sr.register(ref_binary, mov_binary)   # creating transformation matrix \n",
    "out_binary = sr.transform(mov_binary) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "tmat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1,2,figsize = (20,60))\n",
    "axs = axs.ravel()\n",
    "\n",
    "im_reg = np.zeros((2048, 2048,3)) # empty color image\n",
    "im_reg[...,0] = ref_binary\n",
    "im_reg[...,1] = out_binary\n",
    "    \n",
    "im_orig = np.zeros((2048, 2048,3))\n",
    "im_orig[...,0] = ref_binary\n",
    "im_orig[...,1] = mov_binary\n",
    "    \n",
    "axs[0].imshow(im_orig[0:1000, 0:1000])   #before reg\n",
    "    \n",
    "axs[1].imshow(im_reg[0:1000, 0:1000])   #after reg\n",
    "\n",
    "\n",
    "axs[0].title.set_text('Before Registration')\n",
    "axs[1].title.set_text('After Registration')\n",
    "\n",
    "for ax in axs.flat:\n",
    "    ax.set(xlabel='0:1000', ylabel='0:1000')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "reg = np.zeros(mov.shape, dtype=np.uint16) # initialize with the right dtype\n",
    "for Z in range(mov.shape[0]): # Z \n",
    "    print(f\"Z-plane: {Z} registering\")\n",
    "    for ch in range(mov.shape[1]): # channels\n",
    "        reg[Z,ch,...] = sr.transform(mov[Z,ch,...], tmat=tmat)\n",
    "        reg[Z,ch,...] = to_uint16(reg[Z,ch,...])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reg_max = reg.max(0)\n",
    "reg_v = reg_max[0, ...]\n",
    "reg_v = reg_v[100:200, 100:200]\n",
    "reg_v.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# QC by eye, every single Cycle -- pick a random FOV\n",
    "\n",
    "orig = imread('tif/Cycle_9/Cycle_9_F000.tif')\n",
    "orig = orig.max(0)\n",
    "orig = orig[0, ...]\n",
    "orig_v = orig[100:200, 100:200]\n",
    "\n",
    "f, ax = plt.subplots(1,2, figsize = (20,60))\n",
    "ax[0].imshow(orig_v)\n",
    "ax[1].imshow(reg_v)\n",
    "ax[0].title.set_text('Before Registration (DNA Channel)')\n",
    "ax[1].title.set_text('After Registration (DNA Channel)')\n",
    "\n",
    "for ax in ax.flat:\n",
    "    ax.set(xlabel='100:1000', ylabel='100:1000')\n",
    "\n",
    "print(orig_v,\" \\n\" , \"\\n\",reg_v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Registration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check again\n",
    "print('CYCLE_NUMS:', CYCLE_NUMS,'\\n', 'NUM_FOVS:',NUM_FOVS)\n",
    "# CYCLE_NUMS = 10\n",
    "# NUM_FOVS = 214"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(CYCLE_NUMS):\n",
    "    if i == 0:\n",
    "        continue\n",
    "    os.makedirs(f'tmat_Cyc_{i}')\n",
    "    os.makedirs(f'reg_bin_Cyc_{i}')\n",
    "    os.makedirs(f'reg_Cyc_{i}')           "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### CHECK IF REFERENCE FOV IS BEING MATCHED TO SAME MOVED FOV, ACROSS ALL CYCLES AND FOVS IN EVERY ITERATION\n",
    "\n",
    "for c in range(CYCLE_NUMS-1):   \n",
    "    refs = iter(sorted(glob.glob('tif/Cycle_0/*'))) # list of cycle 0 .tif \n",
    "    movs = iter(sorted(glob.glob(f'tif/Cycle_{c+1}/*'))) # cycle 1, 2, 3, .tif list --> FOV000, 001, (002 = error) 005 006 \n",
    "    for FOV in range(0, NUM_FOVS): \n",
    "        #sFOV = str(FOV).zfill(NUM_DIGITS_OF_FOVS)\n",
    "        ref_name = next(refs) \n",
    "        mov_name = next(movs)\n",
    "\n",
    "        ref_num = ref_name.split('_F')[1][0:3]\n",
    "        mov_num = mov_name.split('_F')[1][0:3]\n",
    "        if ref_num != mov_num:\n",
    "            print(\"False\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define jaccard\n",
    "def jaccard(img1, img2):\n",
    "    assert img1.dtype == 'bool', 'input must be boolean'\n",
    "    assert img2.dtype == 'bool', 'input must be boolean'\n",
    "    AND = np.sum(img1&img2)\n",
    "    OR = np.sum(img1|img2)\n",
    "    J = AND/OR\n",
    "    return J\n",
    "reg_J = pd.DataFrame()\n",
    "base_J = pd.DataFrame()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for c in range(2):    \n",
    "    refs = iter(sorted(glob.glob('tif/Cycle_0/*'))) # list of cycle 0 .tif \n",
    "    movs = iter(sorted(glob.glob(f'tif/Cycle_{c+1}/*'))) # cycle 1, 2, 3, .tif list --> FOV000, 001, (002 = error) 005 006 \n",
    "    for FOV in range(0, NUM_FOVS): \n",
    "        #sFOV = str(FOV).zfill(NUM_DIGITS_OF_FOVS)\n",
    "        ref_name = next(refs) \n",
    "        ref = imread(ref_name)\n",
    "        ref = ref.astype(np.uint16)\n",
    "        mov_name = next(movs)\n",
    "        mov = imread(mov_name)\n",
    "        mov = mov.astype(np.uint16)\n",
    "        FOV_num = mov_name.split('_F')[1][0:3]\n",
    "        print(f'cycle {c+1} field {FOV_num} ')\n",
    "\n",
    "        ref_max = ref.max(0)\n",
    "        ref_binary = ref_max[0] > threshold_otsu(ref_max[0]) # nuclei channel\n",
    "        mov_max = mov.max(0)\n",
    "        mov_binary = mov_max[0] > threshold_otsu(mov_max[0]) # nuclei channel\n",
    "        print(\"Got threshold\")\n",
    "        sr = StackReg(StackReg.RIGID_BODY)  \n",
    "        tmat = sr.register(ref_binary, mov_binary) \n",
    "        out = sr.transform(mov_binary) \n",
    "        out = pystackreg.util.to_uint16(out) \n",
    "        \n",
    "#         base_J.loc[FOV_num, str(c+1)] = jaccard(ref_binary, mov_binary)\n",
    "#         reg_J.loc[FOV_num, str(c+1)] = jaccard(ref_binary, out.astype('bool'))\n",
    "\n",
    "        # save binary\n",
    "#         fname_to_save = f'reg_bin_Cyc_{c+1}' + f'/Cycle_{c+1}_F{FOV_num}_bin_reg.tif'\n",
    "#         print('Saving Binary Registered Images...', fname_to_save)\n",
    "#         tifffile.imwrite(fname_to_save, out, imagej=True, photometric = 'minisblack',metadata={'axes':'YX'})\n",
    "\n",
    "        # save tmat\n",
    "        print(\"saving tmat\")\n",
    "        np.save(f'tmat_Cyc_{c+1}' + f'/Cycle_{c+1}_F{FOV_num}_tmat.npy', tmat)\n",
    "        \n",
    "        # THE REGISTRATION STEP\n",
    "#         reg = np.zeros(mov.shape, dtype=np.uint16)               # initialize with the right dtype\n",
    "#         for Z in range(mov.shape[0]):\n",
    "#             print(f\"Z-plane: {Z} registering\")\n",
    "#             for ch in range(mov.shape[1]): \n",
    "#                 reg[Z,ch,...] = sr.transform(mov[Z,ch,...], tmat=tmat)\n",
    "#                 reg[Z,ch,...] = to_uint16(reg[Z,ch,...])\n",
    "\n",
    "#         fname_to_save = f'reg_Cyc_{c+1}' + f'/Cycle_{c+1}_F{FOV_num}_reg.tif'\n",
    "#         print('Saving Registered Images...', fname_to_save)\n",
    "#         tifffile.imwrite(fname_to_save, reg, imagej=True,\n",
    "#                          photometric = 'minisblack',metadata={'axes':'ZCYX'})\n",
    "#base_J.to_csv('base_J.csv')\n",
    "#reg_J.to_csv('reg_J.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Looking at Max and Min Shifts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NUM_FOVS = 195\n",
    "CYCLE_NUMS = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fov_set = set()\n",
    "for idx, c in enumerate(range(CYCLE_NUMS-1)):     \n",
    "    tmats = iter(glob.glob(f'tmat_Cyc_{c+1}/*'))\n",
    "    for sFOV in range(0,NUM_FOVS): \n",
    "        tmat_name = next(tmats)\n",
    "        FOV = tmat_name.split('_F')[1][0:3]\n",
    "        tmat_loaded = np.load(tmat_name)\n",
    "        moveX = tmat_loaded[0,2]\n",
    "        moveY = tmat_loaded[1,2]\n",
    "        if (moveX > 16) | (moveY > 16)|(moveX < -16) | (moveY < -16): # 50 pixels is max\n",
    "            print(c+1, tmat_name, moveX, moveY)\n",
    "            fov_set.add(FOV)\n",
    "print(fov_set)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### QC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# for c in range(CYCLE_NUMS-1):\n",
    "#     refs = iter(sorted(glob.glob('tif/Cycle_0/*'))) # list of cycle 0 .tif \n",
    "#     movs = iter(sorted(glob.glob(f'tif/Cycle_{c+1}/*'))) # cycle 1, 2, 3, .tif list --> FOV000, 001, (002 = error) 005 006 \n",
    "#     for FOV in range(0, NUM_FOVS): # \n",
    "#         #sFOV = str(FOV).zfill(NUM_DIGITS_OF_FOVS)\n",
    "#         ref_name = next(refs) \n",
    "#         ref = imread(ref_name)\n",
    "#         ref = ref.astype(np.uint16)\n",
    "#         mov_name = next(movs)\n",
    "#         mov = imread(mov_name)\n",
    "#         mov = mov.astype(np.uint16)\n",
    "#         FOV_num = mov_name.split('_F')[1][0:3]\n",
    "#         print(f'cycle {c+1} field {FOV_num} ')\n",
    "\n",
    "#         ref_max = ref.max(0)\n",
    "#         ref_binary = ref_max[0] > threshold_otsu(ref_max[0]) # nuclei channel\n",
    "#         mov_max = mov.max(0)\n",
    "#         mov_binary = mov_max[0] > threshold_otsu(mov_max[0]) # nuclei channel\n",
    "#         print(\"Got threshold\")\n",
    "#         sr = StackReg(StackReg.RIGID_BODY)  \n",
    "#         tmat = sr.register(ref_binary, mov_binary) \n",
    "#         out = sr.transform(mov_binary) \n",
    "#         out = pystackreg.util.to_uint16(out) \n",
    "\n",
    "#         base_J.loc[FOV_num, str(c+1)] = jaccard(ref_binary, mov_binary)\n",
    "#         reg_J.loc[FOV_num, str(c+1)] = jaccard(ref_binary, out.astype('bool'))\n",
    "\n",
    "# base_J.to_csv('base_J.csv')\n",
    "# reg_J.to_csv('reg_J.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base_J = pd.read_csv('base_J.csv')\n",
    "base_J = base_J.drop([\"Unnamed: 0\"], axis=1)\n",
    "base_J"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reg_J = pd.read_csv(\"reg_J.csv\")\n",
    "reg_J = reg_J.drop([\"Unnamed: 0\"], axis=1)\n",
    "reg_J"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Making the table for the X and Y shift heatmap   ------ THIS TAKES A VERY SHORT AMOUNT OF TIME\n",
    "\n",
    "for idx, c in enumerate(range(CYCLE_NUMS-1)):\n",
    "    fig, axes = plt.subplots()\n",
    "    tmats = iter(glob.glob(f'tmat_Cyc_{c+1}/*'))\n",
    "    dfX = pd.DataFrame(0, index=['0','1','2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12', '13', '14'], \n",
    "                   columns=['0','1','2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12', '13', '14'])\n",
    "    dfY = pd.DataFrame(0, index=['0','1','2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12', '13', '14'], \n",
    "                   columns=['0','1','2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12', '13', '14'])\n",
    "    for sFOV in range(0,NUM_FOVS):\n",
    "        tmat_name = next(tmats)\n",
    "        tmat_loaded = np.load(tmat_name)\n",
    "        moveX = tmat_loaded[0,2]\n",
    "        moveY = tmat_loaded[1,2]\n",
    "        \n",
    "        col = str(int(tmat_name.split('_')[6]))\n",
    "        row = str(int(tmat_name.split('_')[7][0:3]))\n",
    "\n",
    "        dfX.loc[row, col] = moveX\n",
    "        dfY.loc[row, col] = moveY\n",
    "        \n",
    "    sns.heatmap(dfX, vmin=-20, vmax=20)\n",
    "    pl.suptitle(f\"Cycle {c+1} X-shift\")\n",
    "    \n",
    "#     fig, axes = plt.subplots()\n",
    "#     sns.heatmap(dfY, vmin=-100, vmax=100)\n",
    "#     pl.suptitle(f\"Cycle {c+1} Y-shift\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Making the table for the X and Y shift heatmap   ------ THIS TAKES A VERY SHORT AMOUNT OF TIME\n",
    "\n",
    "for idx, c in enumerate(range(CYCLE_NUMS-1)):\n",
    "    fig, axes = plt.subplots()\n",
    "    tmats = iter(glob.glob(f'tmat_Cyc_{c+1}/*'))\n",
    "    dfX = pd.DataFrame(0, index=['0','1','2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12', '13', '14'], \n",
    "                   columns=['0','1','2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12', '13', '14'])\n",
    "    dfY = pd.DataFrame(0, index=['0','1','2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12', '13', '14'], \n",
    "                   columns=['0','1','2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12', '13', '14'])\n",
    "    for sFOV in range(0,NUM_FOVS):\n",
    "        tmat_name = next(tmats)\n",
    "        tmat_loaded = np.load(tmat_name)\n",
    "        moveX = tmat_loaded[0,2]\n",
    "        moveY = tmat_loaded[1,2]\n",
    "        \n",
    "        col = str(int(tmat_name.split('_')[6]))\n",
    "        row = str(int(tmat_name.split('_')[7][0:3]))\n",
    "\n",
    "        dfX.loc[row, col] = moveX\n",
    "        dfY.loc[row, col] = moveY\n",
    "#     fig, axes = plt.subplots()\n",
    "#     sns.heatmap(dfX, vmin=-20, vmax=20)\n",
    "#     pl.suptitle(f\"Cycle {c+1} X-shift\")\n",
    " \n",
    "    sns.heatmap(dfY, vmin=-20, vmax=20)\n",
    "    pl.suptitle(f\"Cycle {c+1} Y-shift\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfX_shift = pd.DataFrame()\n",
    "dfY_shift = pd.DataFrame()\n",
    "for idx, c in enumerate(range(CYCLE_NUMS-1)):\n",
    "    tmats = iter(sorted(glob.glob(f'tmat_Cyc_{c+1}/*')))\n",
    "    for sFOV in range(0,NUM_FOVS):\n",
    "        tmat_name = next(tmats)\n",
    "        FOV_num = tmat_name.split('_F')[1][0:3]\n",
    "        tmat_loaded = np.load(tmat_name)\n",
    "        moveX = tmat_loaded[0,2]\n",
    "        moveY = tmat_loaded[1,2]\n",
    "\n",
    "        dfX_shift.loc[FOV_num, str(c+1)] = moveX\n",
    "        dfY_shift.loc[FOV_num, str(c+1)] = moveY\n",
    "dfX_shift.to_csv('X_shift.csv')\n",
    "dfY_shift.to_csv('Y_shift.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfY_shift = pd.read_csv(\"Y_shift.csv\")\n",
    "dfY_shift = dfY_shift.drop([\"Unnamed: 0\"], axis=1)\n",
    "dfY_shift\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfX_shift = pd.read_csv(\"X_shift.csv\")\n",
    "dfX_shift = dfX_shift.drop([\"Unnamed: 0\"], axis=1)\n",
    "dfX_shift\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(CYCLE_NUMS-1, 3, figsize=(15,30))\n",
    "\n",
    "for cycle in range(CYCLE_NUMS-1):\n",
    "    axs[cycle][0].hist(reg_J[f'{cycle+1}'], bins = 50, alpha=0.5, label='reg_J', range=[0,1])\n",
    "    axs[cycle][0].hist(base_J[f'{cycle+1}'], bins = 50, alpha=0.5, label='base_J', range=[0,1])\n",
    "    \n",
    "    axs[cycle][1].hist(dfY_shift[f'{cycle+1}'], bins = 50, alpha=0.5, label='dfX_shift', range=[-20,20])\n",
    "    axs[cycle][2].hist(dfX_shift[f'{cycle+1}'], bins = 50, alpha=0.5, label='dfX_shift', range=[-20,20])\n",
    "    \n",
    "    axs[cycle][0].title.set_text(f'Change in Jaccard Index - Cycle {cycle+1}')\n",
    "    axs[cycle][1].title.set_text(f'Y Shift - Cycle {cycle+1}')\n",
    "    axs[cycle][2].title.set_text(f'X Shift - Cycle {cycle+1}')\n",
    "\n",
    "#     for ax in axs.flat:\n",
    "#         ax.set(xlabel='', ylabel='')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# We are trying to look for registered binaries that have intensity percentage less than 0.1 percent (< 0.001)\n",
    "\n",
    "hist = pd.DataFrame()\n",
    "err_lst = []\n",
    "err_lst2 = []\n",
    "err_lst3 = []\n",
    "\n",
    "for c in tqdm(range(CYCLE_NUMS-1)):    \n",
    "    binas = iter(glob.glob(f'reg_bin_Cyc_{c+1}/*')) \n",
    "    for FOV in range(0, NUM_FOVS): # \n",
    "        #sFOV = str(FOV).zfill(NUM_DIGITS_OF_FOVS)\n",
    "        bina_name = next(binas)\n",
    "        bina = imread(bina_name)\n",
    "        bina = bina.astype(np.uint16)\n",
    "        FOV_num = bina_name[-15:-12]\n",
    "        \n",
    "        percentage = (np.sum(bina))/(2048*2048)    \n",
    "        \n",
    "        if percentage < 0.001:\n",
    "            err_lst.append(bina_name)\n",
    "        if percentage < 0.01:\n",
    "            err_lst2.append(bina_name)\n",
    "        if percentage > 0.75:\n",
    "            err_lst3.append(bina_name)\n",
    "            \n",
    "        hist.loc[FOV_num, f'{c+1}'] = percentage    # this is the table of percent of signal in images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "err_lst"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "err_lst\n",
    "uniq_fov = set()\n",
    "for i in err_lst:\n",
    "    FOV = i.split('_F')[1][0:3]\n",
    "    uniq_fov.add(FOV)\n",
    "uniq_fov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "err_lst2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "err_lst2\n",
    "uniq_fov2 = set()\n",
    "for i in err_lst2:\n",
    "    FOV = i.split('_F')[1][0:3]\n",
    "    uniq_fov2.add(FOV)\n",
    "uniq_fov2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "err_lst3\n",
    "uniq_fov3 = set()\n",
    "for i in err_lst3:\n",
    "    FOV = i.split('_F')[1][0:3]\n",
    "    uniq_fov3.add(FOV)\n",
    "uniq_fov3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import pylab as pl\n",
    "for c in range(CYCLE_NUMS-1):\n",
    "    hist.hist(column=f'{c+1}', bins = 80, range=[0, 0.3])\n",
    "    pl.suptitle(f\"Cycle {c+1}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Merging and Cropping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "NUM_FOVS = 195"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## print('X shift and Y shift max min values')\n",
    "x_max_list=[]\n",
    "x_min_list=[]\n",
    "y_max_list=[]\n",
    "y_min_list=[]\n",
    "X_indices =[]\n",
    "Y_indices =[]\n",
    "X_SHIFT_df = pd.DataFrame\n",
    "Y_SHIFT_df = pd.DataFrame\n",
    "\n",
    "for c in range(CYCLE_NUMS-1):     \n",
    "    tmats = iter(glob.glob(f'tmat_Cyc_{c+1}/*'))\n",
    "    X_SHIFT = []\n",
    "    Y_SHIFT = []\n",
    "    X_RIG = []\n",
    "    Y_RIG = []\n",
    "    for sFOV in range(0,NUM_FOVS): \n",
    "        tmat_name = next(tmats)\n",
    "        fov = tmat_name.split('_F')[1][0:3]\n",
    "        tmat_loaded = np.load(tmat_name)\n",
    "        moveX = tmat_loaded[0,2]\n",
    "        moveY = tmat_loaded[1,2]\n",
    "        if moveX > 0:\n",
    "            moveX = moveX + 2048*np.tan(np.arcsin(tmat_loaded[0,1]))\n",
    "            X_RIG = []\n",
    "        if moveY < 0:\n",
    "            moveY = moveY - 2048*np.tan(np.arcsin(tmat_loaded[0,1]))\n",
    "        X_SHIFT.append(moveX)\n",
    "        Y_SHIFT.append(moveY)\n",
    "        \n",
    "#         X_SHIFT_df.loc[fov, c] = moveX\n",
    "#         Y_SHIFT_df.loc[fov, c] = moveY\n",
    "\n",
    "    X_max = max(X_SHIFT)\n",
    "    X_min = min(X_SHIFT)\n",
    "    Y_max = max(Y_SHIFT)\n",
    "    Y_min = min(Y_SHIFT)\n",
    "\n",
    "    print('\\n', f'tmat_Cyc_{c+1} \\n', X_max,X_min,Y_max,Y_min)\n",
    "\n",
    "    def round_shift(val):\n",
    "        import math\n",
    "        if val <0: # if negative,\n",
    "            val = math.floor(val)\n",
    "        else:\n",
    "            val = math.ceil(val)\n",
    "        return val\n",
    "\n",
    "    X_max = round_shift(X_max)\n",
    "    X_min = round_shift(X_min)\n",
    "    Y_max = round_shift(Y_max)\n",
    "    Y_min = round_shift(Y_min)\n",
    "    \n",
    "    x_max_list.append(X_max)\n",
    "    x_min_list.append(X_min)\n",
    "    y_max_list.append(Y_max)\n",
    "    y_min_list.append(Y_min)\n",
    "    print(f'X_max:{X_max}',f'X_min:{X_min}',f'Y_max:{Y_max}',f'Y_min:{Y_min}')\n",
    "\n",
    "X_max_total = max(x_max_list)\n",
    "X_min_total = min(x_min_list)\n",
    "Y_max_total = max(y_max_list)\n",
    "Y_min_total = min(y_min_list)\n",
    "\n",
    "print('\\n', f'X_max_total:{X_max_total}', f'X_min_total:{X_min_total}', \n",
    "      f'Y_max_total:{Y_max_total}', f'Y_min_total:{Y_min_total}')\n",
    "\n",
    "# X_SHIFT_df.to_csv('X_SHIFT_df.csv',sep=',')\n",
    "# Y_SHIFT_df.to_csv('Y_SHIFT_df.csv',sep=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_total = len(mov[0][0][0])    # 2048\n",
    "X_total = len(mov[0][0][1])    # 2048\n",
    "print(Y_total, X_total)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###############################\n",
    "X_abs_shift = [abs(X_max_total), abs(X_min_total)]\n",
    "Y_abs_shift = [abs(Y_max_total), abs(Y_min_total)] \n",
    "\n",
    "def crop(X_x, X_n, Y_x, Y_n, img, pad):\n",
    "    Y_total = img.shape[3]    # 2048\n",
    "    X_total = img.shape[-1]     # 2048\n",
    "    print('image shape', Y_total, X_total)\n",
    "    \n",
    "    if X_x * X_n > 0:   # same signs --> both negative or both positive\n",
    "        Xlength = X_total - (max(X_abs_shift))\n",
    "        print(Xlength)\n",
    "        if X_x > 0: # moving left\n",
    "            img = img[...,:, 0+pad:Xlength-pad]\n",
    "        if X_x < 0:  # moving right \n",
    "            img = img[...,:, abs(X_n)+pad:X_total-pad]\n",
    "\n",
    "    if X_x * X_n < 0:   # diff signs --> one is positive and other is negative\n",
    "        Xlength = X_total  - (abs(X_x) + abs(X_n))\n",
    "        print(Xlength)\n",
    "        img = img[...,:, abs(X_n)+pad: X_total - abs(X_x)-pad] \n",
    "        \n",
    "    if X_x == 0 and X_n < 0:\n",
    "        Xlength = X_total  - abs(X_n)\n",
    "        print(Xlength)\n",
    "        img = img[...,:, abs(X_n)+pad:X_total-pad]\n",
    "    if X_x == 0 and X_n == 0:\n",
    "        Xlength = X_total\n",
    "        print(Xlength)\n",
    "        img = img[...,:, 0+pad:X_total-pad]\n",
    "    if X_x > 0 and X_n == 0:\n",
    "        Xlength = X_total - X_x\n",
    "        print(Xlength)\n",
    "        img = img[...,:, 0+pad:Xlength-pad]\n",
    "### Y ## #\n",
    "    if Y_x * Y_n > 0:   \n",
    "        Ylength = Y_total - (max(Y_abs_shift))\n",
    "        print(Ylength)\n",
    "        if Y_x > 0: # up\n",
    "            img = img[..., 0+pad:Ylength-pad, :]\n",
    "        if Y_x < 0: # down\n",
    "            img = img[...,abs(Y_n)+pad:Y_total-pad, :]\n",
    "            \n",
    "    if Y_x * Y_n < 0:  # Y_x > 0 , Y_n < 0\n",
    "        Ylength = Y_total  - (abs(Y_x) + abs(Y_n)) # \n",
    "        print(Ylength)\n",
    "        img = img[...,abs(Y_n)+pad:Y_total-abs(Y_x)-pad, :] \n",
    "        \n",
    "    if Y_x == 0 and Y_n < 0: # down\n",
    "        Ylength = Y_total  - abs(Y_n)\n",
    "        print(Ylength)\n",
    "        img = img[..., abs(Y_n)+pad:Y_total-pad, :]\n",
    "        \n",
    "    if Y_x == 0 and Y_n == 0:\n",
    "        Ylength = Y_total\n",
    "        print(Ylength)\n",
    "        img = img[..., 0+pad:Y_total-pad,:]\n",
    "        \n",
    "    if Y_x > 0 and Y_n == 0:\n",
    "        Ylength = Y_total - Y_x ## \n",
    "        print(Ylength)\n",
    "        img = img[..., 0+pad:Ylength-pad, :] ## up  \n",
    "    Xlength = Xlength - (pad*2)\n",
    "    Ylength = Ylength - (pad*2)\n",
    "    return Xlength, Ylength, img"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "IN_DIR = 'reg'\n",
    "REF_DIR = 'tif' # cycle0\n",
    "MERGE_DIR = 'merged' # output directory to save\n",
    "Z = 13\n",
    "final_ch = 38\n",
    "\n",
    "pad = 10 # 5 pixel padding\n",
    "\n",
    "!mkdir merged"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "refs = iter(glob.glob('tif/Cycle_0/*')) \n",
    "for FOV in range(NUM_FOVS):\n",
    "    ref_name = next(refs)\n",
    "    FOV_num = ref_name.split('_F')[1][0:3]\n",
    "    print(\"FOV\", FOV_num)\n",
    "    #sFOV = str(FOV).zfill(NUM_DIGITS_OF_FOVS)\n",
    "    img = imread(f'{REF_DIR}/Cycle_0/Cycle_F{FOV_num}.tif') # reference image that did not move\n",
    "    img = img.astype(np.uint16)\n",
    "    print(\"image max:\",img.max())\n",
    "    for cycle in range(CYCLE_NUMS-1):    \n",
    "        fname = f'{IN_DIR}_Cyc_{cycle+1}/Cycle_{cycle+1}_F{FOV_num}_{IN_DIR}.tif'\n",
    "        print('Appending... ', fname)\n",
    "        im_to_add = imread(fname).astype(np.uint16)\n",
    "        print(im_to_add.max())\n",
    "        im_to_add = im_to_add[:,:,...] \n",
    "        img = np.append(img, im_to_add, axis=1) # concatenate along channel index\n",
    "        print(f\"Shape = {img.shape}\") \n",
    "    fname = f'/F{FOV_num}.tif'\n",
    "    print('saving', './{MERGE_DIR}'+fname)\n",
    "    print('dtype of ', img.dtype)\n",
    "\n",
    "    ########## CROP ################ --> CHANGE everytime depending on shifts\n",
    "    Xlength, Ylength, img = crop(X_max_total, X_min_total, Y_max_total, Y_min_total, img, pad)\n",
    "\n",
    "    assert img.shape[3] == Xlength, \"Check X size\"\n",
    "    assert img.shape[2] == Ylength, \"Check Y size\"\n",
    "\n",
    "    ### FINAL CHECK before saving ### \n",
    "    assert img.shape[0] == Z, \"Check ZCYX\"\n",
    "    assert img.shape[1] == final_ch, \"check final merge size\"\n",
    "    print(\"image max:\", img.max())\n",
    "\n",
    "    tifffile.imwrite(\n",
    "        f'./{MERGE_DIR}'+fname,\n",
    "        img,\n",
    "        imagej=True,\n",
    "        photometric='minisblack',\n",
    "        metadata={'axes': 'ZCYX'},\n",
    "    )\n",
    "\n",
    "    del img # clear memory\n",
    "    del im_to_add\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NUM_FOVS = 195"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "int_counts = pd.DataFrame()\n",
    "one_to_five_list = []\n",
    "sixfivek_list = []\n",
    "z_list = []\n",
    "\n",
    "merged = iter(glob.glob('merged/*')) \n",
    "for FOV in range(NUM_FOVS):\n",
    "    merged_name = next(merged)\n",
    "    img = imread(merged_name)\n",
    "    img = img.astype(np.uint16)\n",
    "    FOV_num = merged_name.split('/F')[1][0:3]\n",
    "    print(\"FOV\", FOV_num)\n",
    "    \n",
    "    overall_count_zero = np.count_nonzero(img == 0)\n",
    "    x = np.count_nonzero((0 < img) & (img < 6))\n",
    "    y = np.count_nonzero(65000 < img)\n",
    "     \n",
    "    for Z in range(img.shape[0]):\n",
    "        print(Z)\n",
    "        for ch in range(img.shape[1]): \n",
    "            count_zeros = np.count_nonzero(img[Z,ch,...] == 0)\n",
    "            if count_zeros > 0:\n",
    "                z_list.append((FOV_num, Z, ch))\n",
    "            count = np.count_nonzero((0 < img[Z,ch,...]) & (img[Z,ch,...] < 6))\n",
    "            if count > 0:\n",
    "                one_to_five_list.append((FOV_num, Z, ch))\n",
    "            count = np.count_nonzero(65000 < img[Z,ch,...])\n",
    "            if count > 0:\n",
    "                sixfivek_list.append((FOV_num, Z, ch))\n",
    "                \n",
    "    int_counts.loc[FOV_num, 'FOV_num'] = FOV_num\n",
    "    int_counts.loc[FOV_num, 'total Zero count'] = overall_count_zero\n",
    "    int_counts.loc[FOV_num, 'one to five'] = x\n",
    "    int_counts.loc[FOV_num, 'greater 65K'] = y\n",
    "int_counts.to_csv('pixel_intensity_counts_table.csv')    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "int_counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(sixfivek_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(one_to_five_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(z_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for idx, val in enumerate(int_counts['one to five']):\n",
    "    if val != 0:\n",
    "        print(int(int_counts.iloc[idx, 2]), \"FOV:\", int_counts.iloc[idx, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_FOV_set = set()\n",
    "for i in z_list:\n",
    "    z_FOV_set.add(i[0])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# LOOKING FOR WHICH FOVS HAVE THE BLACK BOXES\n",
    "black_box = []\n",
    "for idx, val in enumerate(int_counts['total Zero count']):\n",
    "    if val > 2000:\n",
    "        black_box.append(int_counts.iloc[idx, 0])\n",
    "        print(int_counts.iloc[idx, 0], val)\n",
    "black_box"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "len(black_box)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "one_to_five_set = set()\n",
    "for i in one_to_five_list:\n",
    "    one_to_five_set.add(i[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "len(one_to_five_set)               # num FOVs that have pixel intensity 1-5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(one_to_five_list)    # number of single images that contains -- surrounding black box, black circle artifacts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "one_to_five_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Checking Location of Pixels in Images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# 0 pixel image \n",
    "for i in z_list:            \n",
    "    if i[0] == '122':\n",
    "        print(f\"Z:{i[1]}\", \" \",f\"ch:{i[2]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = 'merged/F000.tif'\n",
    "img = imread(img)\n",
    "np.where(img == 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# >65k pixel image\n",
    "for i in sixfivek_list:           \n",
    "    if i[0] == '002':\n",
    "        print(f\"Z:{i[1]}\", \" \",f\"ch:{i[2]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "img = 'merged/F002.tif'\n",
    "img = imread(img)\n",
    "np.where(65000 < img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1-5 pixel image\n",
    "for i in one_to_five_list:            \n",
    "    if i[0] == '000':\n",
    "        print(f\"Z:{i[1]}\", \" \",f\"ch:{i[2]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = 'merged/F000.tif'\n",
    "img = imread(img)\n",
    "np.where((0 < img) & (img < 6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df_black_box = pd.DataFrame()\n",
    "for idx, i in enumerate(z_list):\n",
    "    if i[0] in black_box:\n",
    "        df_black_box.loc[idx, 'FOV_num'] = i[0]\n",
    "        df_black_box.loc[idx, 'Z'] = i[1]\n",
    "        df_black_box.loc[idx, 'channel'] = i[2]\n",
    "        \n",
    "df_black_box\n",
    "df_black_box.to_csv('black_box_locations.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# total num of pixels is 3995601 (because 2019*1979)\n",
    "total = img.shape[2] * img.shape[3]\n",
    "# here we get percent of 0 intensity count, 0<x<6 and greater than 65k\n",
    "percents = int_counts.copy()\n",
    "percents['total Zero count'] = percents['total Zero count'].div(total)\n",
    "percents['one to five'] = percents['one to five'].div(total)\n",
    "percents['greater 65K'] = percents['greater 65K'].div(total)\n",
    "percents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "percents.to_csv('percent_intensity_table.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "black_circles = percents.loc[(percents['one to five'] > 0)&(percents['greater 65K'] > 0)&(percents['total Zero count'] < 0.00050055048)]  \n",
    "black_circles\n",
    "# THESE MUST BE THE BLACK CIRCLE ARTIFACT ONES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "black_circles.to_csv('black_circle_pixel_intensity_table.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "black_circle_lst = []\n",
    "for i in black_circles['FOV_num']:\n",
    "    black_circle_lst.append(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "black_circle_lst"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_black_circle = pd.DataFrame()\n",
    "count = 0\n",
    "for idx, i in enumerate(one_to_five_list):\n",
    "    if i[0] in black_circle_lst:\n",
    "        df_black_circle.loc[idx, 'FOV_num'] = i[0]\n",
    "        df_black_circle.loc[idx, 'Z'] = i[1]\n",
    "        df_black_circle.loc[idx, 'channel'] = i[2]\n",
    "        count += 1\n",
    "        \n",
    "df_black_circle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_black_circle.to_csv('black_circle_locations.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = imread('merged/F000.tif')\n",
    "image = image[0, 32, ...]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "t = 0\n",
    "binary_mask = image <= t\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(7, 7), sharex=True, sharey=True)\n",
    "plt.title('Zero Intensity Masked Image')\n",
    "plt.imshow(binary_mask, cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# BINARY MASKING FOR 0 INTENSITY PIXELS -- ANOTHER WAY OF FINDING LOCATION OF BLACK BOXES\n",
    "validation = []\n",
    "t = 0\n",
    "merged = iter(glob.glob('merged/*')) \n",
    "for FOV in tqdm(range(NUM_FOVS)):\n",
    "    merged_name = next(merged)\n",
    "    img = imread(merged_name)\n",
    "    img = img.astype(np.uint16)\n",
    "    FOV_num = merged_name.split('/F')[1][0:3]\n",
    "    print(\"FOV\", FOV_num)\n",
    "\n",
    "    for Z in range(img.shape[0]):\n",
    "        for ch in range(img.shape[1]): \n",
    "            binary_mask = img[Z,ch,...] <= t\n",
    "            if binary_mask.sum() > 2000:\n",
    "                validation.append((FOV_num, Z, ch))\n",
    "                \n",
    "done = open(\"done.txt\", \"a\")\n",
    "done.write(\"done\")\n",
    "done.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df_validation = pd.DataFrame(validation)\n",
    "df_validation = df_validation.rename(columns={0: \"FOV_num\", 1: \"Z\", 2:'channel'})\n",
    "df_validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_validation.to_csv('black_box_locations_correct.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_validation = pd.read_csv('black_box_locations_correct.csv')\n",
    "df_validation = df_validation.drop(['Unnamed: 0'], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "validation_lst = list(df_validation.itertuples(index=False, name=None))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for idx, val in tqdm(enumerate(validation_lst)):\n",
    "    FOV = str(val[0]).zfill(3)\n",
    "    img = imread(f'merged/F{FOV}.tif')\n",
    "    Z = val[1]\n",
    "    ch = val[2]\n",
    "    \n",
    "    cnt = np.count_nonzero(img[Z, ch, ...] == 0)\n",
    "    ttl = img.shape[2]*img.shape[3]\n",
    "    df_validation.loc[idx, 'percent_bbox'] = (cnt/ttl)*100\n",
    "df_validation.to_csv('bbox_locs_and_percs.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df_validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "validation_FOVs = set()\n",
    "for i in validation:\n",
    "    validation_FOVs.add(i[0])\n",
    "validation_FOVs "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "validation_FOVs = np.array(sorted(list(validation_FOVs)))\n",
    "black_box_arr = np.array(black_box)\n",
    "np.array_equal(validation_FOVs, black_box_arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = imread('merged/F000.tif')\n",
    "arr = arr[0, 32, ...]\n",
    "#plt.imshow(arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len((np.where(arr == 0))[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "percent = 0\n",
    "while np.percentile(arr, percent) == 0:\n",
    "    print(np.percentile(arr, percent), percent)\n",
    "    percent += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NUM_FOVS = 211"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "percentiles = pd.DataFrame()\n",
    "\n",
    "merged = iter(glob.glob('merged/*')) \n",
    "for FOV in tqdm(range(NUM_FOVS)):\n",
    "    merged_name = next(merged)\n",
    "    img = imread(merged_name)\n",
    "    img = img.astype(np.uint16)\n",
    "    FOV_num = merged_name.split('/F')[1][0:3]\n",
    "    print(\"FOV\", FOV_num)\n",
    "\n",
    "    for ch in range(img.shape[1]): \n",
    "        arr = img[:,ch,...]\n",
    "        percent = 0\n",
    "        while np.percentile(arr, percent) == 0:\n",
    "            percent += 1\n",
    "        percentiles.loc[FOV_num, f'ch{ch}_cutoff_perc'] = percent # by the +1, you will know if there are ANY 0's\n",
    "        \n",
    "        percentiles.loc[FOV_num, f'ch{ch}_0th_perc'] = np.percentile(arr, 0)\n",
    "        percentiles.loc[FOV_num, f'ch{ch}_0.001st_perc'] = np.percentile(arr, 0.001)\n",
    "        percentiles.loc[FOV_num, f'ch{ch}_0.01st_perc'] = np.percentile(arr, 0.01)\n",
    "        percentiles.loc[FOV_num, f'ch{ch}_0.1st_perc'] = np.percentile(arr, 0.1)\n",
    "        percentiles.loc[FOV_num, f'ch{ch}_0.5th_perc'] = np.percentile(arr, 0.5)\n",
    "        percentiles.loc[FOV_num, f'ch{ch}_1st_perc'] = np.percentile(arr, 1)\n",
    "        percentiles.loc[FOV_num, f'ch{ch}_5th_perc'] = np.percentile(arr, 5)\n",
    "        percentiles.loc[FOV_num, f'ch{ch}_10th_perc'] = np.percentile(arr, 10)\n",
    "        percentiles.loc[FOV_num, f'ch{ch}_90th_perc'] = np.percentile(arr, 90)\n",
    "        percentiles.loc[FOV_num, f'ch{ch}_95th_perc'] = np.percentile(arr, 95)\n",
    "        percentiles.loc[FOV_num, f'ch{ch}_99th_perc'] = np.percentile(arr, 99)\n",
    "        percentiles.loc[FOV_num, f'ch{ch}_100th_perc'] = np.percentile(arr, 100)\n",
    "\n",
    "\n",
    "percentiles.to_csv('percentiles.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "percentiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "zeroth_perc = percentiles.filter(regex=(\".*_0th.*\"))\n",
    "zeroth_perc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# histogram of maz 0th percentile values across all channels\n",
    "\n",
    "plt.hist(zeroth_perc.max().to_frame()[0],bins = 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Channel 35\n",
    "plt.hist(percentiles['ch35_0th_perc'], bins = 100, alpha=0.5, range=[0, 420], label = '0_perc')\n",
    "plt.hist(percentiles['ch35_0.001st_perc'], bins = 100, alpha=0.5, range=[0, 420], label = '0.001_perc')\n",
    "plt.hist(percentiles['ch35_0.1st_perc'], bins = 100, alpha=0.5, range=[0, 420], label = '0.1_perc')\n",
    "plt.hist(percentiles['ch35_5th_perc'], bins = 100, alpha=0.5, range=[0, 420], label = '5_perc')\n",
    "plt.hist(percentiles['ch35_10th_perc'], bins = 100, alpha=0.5, range=[0, 420], label = '10_perc')\n",
    "plt.legend(loc='upper left')\n",
    "plt.suptitle(f\"Channel 35 -- percentiles\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Channel 36\n",
    "plt.hist(percentiles['ch36_0th_perc'], bins = 100, alpha=0.5, range=[0, 420], label = '0_perc')\n",
    "plt.hist(percentiles['ch36_0.001st_perc'], bins = 100, alpha=0.5, range=[0, 420], label = '0.001_perc')\n",
    "plt.hist(percentiles['ch36_0.1st_perc'], bins = 100, alpha=0.5, range=[0, 420], label = '0.1_perc')\n",
    "plt.hist(percentiles['ch36_5th_perc'], bins = 100, alpha=0.5, range=[0, 420], label = '5_perc')\n",
    "plt.hist(percentiles['ch36_10th_perc'], bins = 100, alpha=0.5, range=[0, 420], label = '10_perc')\n",
    "plt.legend(loc='upper left')\n",
    "plt.suptitle(f\"Channel 36 -- percentiles\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for i in tqdm(range(38)):\n",
    "    fig, axes = plt.subplots()\n",
    "    plt.hist(percentiles[f'ch{i}_0th_perc'], bins = 100, alpha=0.5, range=[0, 110], label = '0_perc')\n",
    "    plt.hist(percentiles[f'ch{i}_0.001st_perc'], bins = 100, alpha=0.5, range=[0, 110], label = '0.001_perc')\n",
    "    plt.hist(percentiles[f'ch{i}_0.1st_perc'], bins = 100, alpha=0.5, range=[0, 110], label = '0.1_perc')\n",
    "    plt.hist(percentiles[f'ch{i}_5th_perc'], bins = 100, alpha=0.5, range=[0, 110], label = '5_perc')\n",
    "    plt.hist(percentiles[f'ch{i}_10th_perc'], bins = 100, alpha=0.5, range=[0, 110], label = '10_perc')\n",
    "    plt.legend(loc='upper left')\n",
    "    plt.suptitle(f\"Channel {i} -- percentiles\")\n",
    "    \n",
    "#     fig, axes = plt.subplots()\n",
    "#     sns.heatmap(dfY, vmin=-100, vmax=100)\n",
    "#     pl.suptitle(f\"Cycle {c+1} Y-shift\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cutoffs = percentiles.filter(regex=(\".*cutoff.*\"))\n",
    "cutoffs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tup = []\n",
    "for index, row in cutoffs.iterrows():\n",
    "    for col, val in enumerate(row):\n",
    "        if val > 0:\n",
    "            tup.append((index, col, val))\n",
    "tup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bb_perc = pd.DataFrame(tup)\n",
    "bb_perc = bb_perc.rename(columns={0: \"FOV_num\", 1: \"channel\", 2:'percentage'})\n",
    "bb_perc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bb_perc.to_csv('bbox_percentages.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bc = pd.read_csv('black_circle_locations.csv')\n",
    "bc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = set()\n",
    "for i in bc['FOV_num']:\n",
    "    s.add(str(i).zfill(3))\n",
    "s = pd.DataFrame(sorted(s))\n",
    "s.to_csv('b_circle_FOVs.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test = imread('merged/F132.tif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.where(test[:, 35, ...] == 0)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(x[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test = imread('merged/F002.tif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.where(test[:, 28, ...] == 0)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(x[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as patches\n",
    "from pathlib import Path\n",
    "import numpy as np\n",
    "import os\n",
    "import sys\n",
    "import glob\n",
    "from imageio import volread as imread\n",
    "\n",
    "from skimage.segmentation import expand_labels\n",
    "import tifffile\n",
    "from skimage.segmentation import *\n",
    "from skimage import measure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mkdir mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAGxCAYAAAD/MbW0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABFAUlEQVR4nO3deVxUVf8H8M/MAIOK4oKAIooL7rmkiWipFIobZampqRipGUmLPL9SXCCXxOp51EyKVyqhloqZmaSZiWuJj7nlkjsgpUKYioYyKHN+f/hyHmdYZi7cYWa4n3ev+3p175n53u8MOPPlnHPPVQkhBIiIiEix1LZOgIiIiGyLxQAREZHCsRggIiJSOBYDRERECsdigIiISOFYDBARESkciwEiIiKFYzFARESkcCwGiIiIFI7FABERkcKxGCAiIrITe/fuRWhoKBo2bAiVSoVNmzaZfc7u3bvx+OOPQ6vVokWLFkhKSpJ8XhYDREREdiI/Px8dO3ZEfHy8RY/PyMjAoEGDEBQUhGPHjuHtt9/GhAkT8OOPP0o6r4o3KiIiIrI/KpUK3377LYYMGVLqY6ZOnYotW7bg5MmThmMjR47EzZs3sW3bNovPxZ4BIiIiK9LpdLh165bRptPpZImdlpaG4OBgo2MhISFIS0uTFMdJlmyIiIiqkHvX0mWLFbd0FWbPnm10LDY2Fu+9916FY2dnZ8PLy8vomJeXF27duoW7d++iWrVqFsWxq2Kgom++s0czAICTi0+F4twvvGyX+VTVOPb2PttbPlU1TlV9n/m6SvbwdcnxJfswlqOIjo5GVFSU0TGtVmujbEpmV8UAERGRXdAXyRZKq9Va7cvf29sbOTk5RsdycnJQq1Yti3sFABYDRERExQm9rTOwSGBgILZu3Wp07KeffkJgYKCkOJxASEREZCf++ecfHDt2DMeOHQPw4NLBY8eOISsrC8CDIYewsDDD41977TWkp6fj3XffxZkzZ/Dpp59i/fr1mDJliqTzsmeAiIjIlN42PQOHDh1CUFCQYf/hXINx48YhKSkJV69eNRQGANC0aVNs2bIFU6ZMwccff4xGjRph+fLlCAkJkXReFgNEREQmhI2GCfr06YOylv8paXXBPn364OjRoxU6L4sBIiIiUzbqGbAVzhkgIiJSOPYMEBERmXKQqwnkwmKAiIjIlIzrDDgCDhMQEREpHHsGiIiITHGYQLr79++joKAAbm5ucoQjIiKyLV5NULqUlJRi1zi+//77cHNzQ+3atdGvXz/cuHFDzvyIiIjIyiQVAwsXLkR+fr5hf//+/YiJicGsWbOwfv16/PHHH5g7d67sSRIREVUmIfSybY5A0jDBqVOnsHDhQsP+hg0b0LdvX8yYMQMA4OrqirfeesvoMURERA6HwwSlu337NurVq2fY//nnn/HMM88Y9tu1a4crV67Ilx0RERFZnaRiwMfHB6dPnwbw4M5Kv/32G3r06GFo//vvv1G9enV5MyQiIqpsQi/f5gAkDRMMHz4cb7/9NqZPn46tW7fC29sb3bt3N7QfOnQIrVq1kj1JIiKiSqWwRYckFQMxMTG4fPky3nzzTXh7e+PLL7+ERqMxtK9duxahoaGyJ0lERFSpHOQverlIKgaqVauGVatWldq+a9euCidERERElUvSnIHY2Fjs3bsXhYWF1sqHiIjI9vR6+TYHIKkYWLVqFfr06YPatWvjmWeewbx58/DLL7/g/v371sqPiIio8ilsAqGkYiAjIwPp6emIj49Ho0aNsHz5cjz11FOoU6cO+vfvjw8++AAHDx60Vq5ERERkBZLvWujn54fw8HCsXLkSmZmZuHjxIj7++GN4enpi/vz5RpcaEhEROSSFDRNU6EZFly5dwt69e7Fnzx7s3bsX9+7dQ69evcw+T6fTQafTGR3TarW8nzIREdkFIZR1aaGk79+srCysWrUK4eHhaNq0Kdq3b481a9agVatW+PLLL3Hz5k3s3LnTbJy4uDi4u7sbbXFxceV+EURERFR+knoG/Pz80LhxY0RERCAiIgJdunQxWmfAUtHR0YiKijI6ptVqgduXJcciIiKSnYNM/JOLpGLgxRdfxJ49e/DBBx/gl19+Qe/evREUFITOnTtDpVJZHEer1T748jdx77aUbIiIiKzEQcb65SKpGFi3bh0A4MyZM9i1axd2796Njz76CAUFBXjyySfRu3dv9OnTB0888YRVkiUiIqoUCusZKNecvdatWyMiIgLJycnIzs7G/v370alTJ8ybNw+BgYFy50hERERWVO6rCXJycrB7927s3r0bu3btwrlz56DVavHUU0/JmR8REVHl442KSrd+/XpDAXD27Fk4OzvjiSeewIsvvoigoCD06NGjxLkAREREDkVhwwSSioExY8aga9eueP755xEUFISePXuiWrVq1sqNiIiIKoGkYuDGjRuoUaOGtXIhIiKyD7yaoHQ1a9Y0ewmhSqXijYuIiMixcZigdBs3biy1GEhLS8OSJUugV1g1RURE5OgkFQNDhgwpduzs2bOYNm0aUlJSMHr0aMyZM0eu3IiIiGxDYX/YlvveQFeuXMHEiRPx2GOP4f79+zh27BhWrlyJJk2ayJkfERFR5VPYXQslFwN5eXmYOnUqWrRogVOnTiE1NRUpKSlo3769NfIjIiIiK5M0TPDhhx/igw8+gLe3N9auXYvnnnvOWnkRERHZjNJuYSypGJg2bRqqVauGFi1aYOXKlVi5cmWJj9u4caMsyREREdmEg3Tvy0VSMRAWFibp7oREREQOiZcWli4pKclKaRAREZGtlPtGRURERFUWhwmIiIgUTmHDBOVeZ4CIiIiqBvYMEBERmeIwARERkcJxmICIiIiUhD0DREREphQ2TKASQghbJ0FERGRP7m5ZLFusaoPeli2WtXCYgIiISOHsapjg3rX0Cj3f2aMZAMDJxadCce4XXmY+FuTDOGXHsbefl73lY2+vq6rGsbefV0XjPBrLqhQ2gdCuigEiIiK7oLA5AywGiIiITCmsZ4BzBoiIiBSOPQNERESmOExARESkcBwmICIiIiVhzwAREZEpDhMQEREpHIuB0t26dQu1atUCAGzduhX37983tGk0GgwaNEje7IiIiMjqLC4Gvv/+e8yaNQtHjx4FAIwYMQL5+fmGdpVKheTkZAwbNkz+LImIiCqTwm7bY/EEws8//xxvvPGG0bELFy5Ar9dDr9cjLi4OiYmJsidIRERU6fR6+TYHYHExcOLECfTs2bPU9gEDBuDQoUOyJEVERESVx+JhgqtXr0Kr1Rr2d+3aBV9fX8O+m5sb8vLy5M2OiIjIFhzkL3q5WNwzULduXVy4cMGw37VrVzg7Oxv2z58/j7p168qbHRERkS0IvXybA7C4GOjVqxeWLFlSavuSJUvQq1cvWZIiIiKyKc4ZKNnUqVOxfft2DB8+HL/++ivy8vKQl5eHgwcPYujQodixYwemTp1qzVyJiIiqvPj4ePj5+cHV1RUBAQE4ePBgmY9fvHgxWrVqhWrVqsHX1xdTpkxBQUGBpHNaPGegc+fOSE5OxoQJE7Bx40ajtjp16mDdunV4/PHHJZ2ciIjILtno0sLk5GRERUUhISEBAQEBWLx4MUJCQnD27Fl4enoWe/yaNWswbdo0JCYmokePHjh37hxefvllqFQqLFy40OLzSlp06LnnnkPfvn3x448/4vz58wAAf39/9OvXDzVq1JASioiIyH7ZqHt/4cKFmDhxIsLDwwEACQkJ2LJlCxITEzFt2rRij9+/fz969uyJl156CQDg5+eHUaNG4b///a+k81pcDCxduhRjx46Fu7s7nn/+eUknISIiUiqdTgedTmd0TKvVGl2hBwCFhYU4fPgwoqOjDcfUajWCg4ORlpZWYuwePXrgyy+/xMGDB9GtWzekp6dj69atGDt2rKQcLZ4zMGPGDDRo0AAvvfQSdu7cKekkREREDkXGCYRxcXFwd3c32uLi4oqd8tq1aygqKoKXl5fRcS8vL2RnZ5eY5ksvvYQ5c+bgySefhLOzM5o3b44+ffpg+vTpkl6uxcVAdnY2EhIScPXqVfTt2xdNmzbF3Llz8ccff0g6IRERkd2T8dLC6Ohow6T7h9ujf/1XxO7duzF//nx8+umnOHLkCDZu3IgtW7Zg7ty5kuJYXAxUq1YNYWFh2LVrF86fP4+xY8dixYoVaNq0Kfr374+vv/4a9+7dk/xCiIiIqjKtVotatWoZbaZDBADg4eEBjUaDnJwco+M5OTnw9vYuMfasWbMwduxYTJgwAY899hief/55zJ8/H3FxcdBLmPdgcTHwqGbNmmHOnDnIyMjADz/8gHr16uHll1+Gj49PecIRERHZFaEXsm2WcnFxQZcuXZCammo4ptfrkZqaisDAwBKfc+fOHajVxl/lGo3mwWuQcEWEpKsJTKlUKjg5OUGlUkEIYXHPQGmTKcpVmRAREcnNRlcTREVFYdy4cejatSu6deuGxYsXIz8/33B1QVhYGHx8fAxzDkJDQ7Fw4UJ07twZAQEBuHDhAmbNmoXQ0FBDUWCJchUDf/zxB7744gskJSUhKysLvXr1wrJlyzB06FCLnh8XF4fZs2cbHYuNjcWMyLDypENERFQljBgxArm5uYiJiUF2djY6deqEbdu2GSYVZmVlGfUEzJw5EyqVCjNnzsTly5dRv359hIaG4v3335d0XouLgcLCQmzcuBGJiYnYuXMnGjRogHHjxuGVV15Bs2bNJJ00OjoaUVFRRse0Wi1w+7KkOERERFZhw3sKREZGIjIyssS23bt3G+07OTkhNjYWsbGxFTqnxcWAt7c37ty5g8GDByMlJQUhISHFxiksVdL1lQBw73a5whEREclLwlh/VWBxMTBz5kyMHTsW9evXt2Y+REREtucgNxiSi8XFgGm3PgAUFBQgOTkZ+fn56Nu3L/z9/WVNjoiIiKxPUjFw7949fPLJJwAezCEIDAzEqVOnUL16dbz77rv46aefSr38gYiIyGEorGfA4kH/7du3o2/fvob9r776CpcuXcL58+dx48YNDB8+HPPmzbNKkkRERJVKCPk2B2BxMZCVlYW2bdsa9rdv345hw4ahSZMmUKlUeOutt3D06FGrJElERETWY3ExoFarjVYzOnDgALp3727Yr127Nm7cuCFvdkRERLYg442KHIHFxUCbNm2QkpICADh16hSysrIQFBRkaL906VKxOy0RERE5JL2Qb3MAFk8gfPfddzFy5Ehs2bIFp06dwsCBA9G0aVND+9atW9GtWzerJElERETWY3Ex4O/vjx9++AEpKSno168f3njjDaP26tWr4/XXX5c9QSIiokpnwxUIbcHiYqBDhw544oknMH78eIwaNQrVq1c3aq/oUohERER2w0G69+Vi8ZyBPXv2oF27dvi///s/w30J9u3bZ83ciIiIqBJYXAw89dRTSExMxNWrV/HJJ58gMzMTvXv3RsuWLfHBBx8gOzvbmnkSERFVGqHXy7Y5Asl3GqpRowbCw8OxZ88enDt3DsOHD0d8fDwaN26MZ5991ho5EhERVS5eTWC5Fi1aYPr06WjSpAmio6OxZcsWufIiIiKyHU4gtMzevXuRmJiIb775Bmq1Gi+++CLGjx8vZ25ERERUCSQVA1euXEFSUhKSkpJw4cIF9OjRA0uWLMGLL76IGjVqWCtHIiKiyuUg3ftysbgYGDBgAHbs2AEPDw+EhYXhlVdeQatWrayZGxERkW04yMQ/uVhcDDg7O2PDhg0YPHgwNBqNNXMiIiKiSmRxMbB582Zr5kFERGQ/OExARESkcAq7mkDyOgNERERUtbBngIiIyBSHCYiIiJTNUZYRlotKCKGs8oeIiMiMf6KHyhbLLe4b2WJZC3sGiIiITHGYwHbuXUuv0POdPZrJGsfJxadCce4XXpY1jr29LuZTdj5VNU5VfZ/5ukpmb/9OH83JqlgMEBERKRwvLSQiIiIlYc8AERGRKQ4TEBERKZtQWDHAYQIiIiKFY88AERGRKYX1DLAYICIiMqWwFQg5TEBERKRw7BkgIiIyxWECIiIihVNYMcBhAiIiIoWTVAxcvHgRr7zyimG/cePGqFu3rmGrX78+zp49K3uSRERElUkIIdvmCCQNE3zyySfw8vIy7N+4cQMxMTHw9PQEACQnJ2PRokVISEiQN0siIqLKpLBhAknFQGpqKlasWGF0bOjQoWjW7MEdpPz8/DBhwgT5siMiIrIFhRUDkoYJMjMz0bBhQ8P+hAkT4O7ubtj38/PDn3/+KV92REREZHWSegbUajWuXLmCRo0aAQAWLVpk1J6TkwNnZ2f5siMiIrIB3pugDO3atcOOHTtKbf/xxx/Rvn37CidFRERkU3oh3+YAJBUD4eHheP/997Fly5ZibSkpKViwYAHCw8NlS46IiIisT9IwwcSJE7Fz506EhoaidevWaNWqFQDg7NmzOHv2LIYOHYqJEydaJVEiIqJKo6xbE0hfdGjt2rVYs2YNWrZsaSgC/P398dVXX2H9+vXWyJGIiKhSCb2QbXME5VqOeOTIkRg5cmSx43q9Hlu3bsXgwYMrnBgRERFVDlnuTXDhwgUkJiYiKSkJubm5uHfvnhxhiYiIbMNB/qKXS7nvTXD37l2sWrUKvXr1QqtWrbB//37ExMRwnQEiInJ8ehk3ByC5Z+DXX3/F8uXLsW7dOjRv3hyjR4/G/v378emnn6Jt27YWxdDpdNDpdEbHtFot75pERERkA5K+fzt06IDhw4ejXr162L9/P44cOYJ//etfUKlUkk4aFxcHd3d3oy0uLk5SDCIiImvhBMIynD17FiNGjEBQUJDFvQAliY6ORlRUlNExrVYL3L5c7phERESycZDufblIKgbS09ORlJSEiIgI3L17F6NGjcLo0aMl9wxotdoHX/4m7t2WFIaIiMgqHOUverlIGibw8fHBjBkzcOHCBaxevRrZ2dno2bMn7t+/j6SkJJw7d85aeRIREZGVlHvO3tNPP40vv/wSV69exdKlS7Fz5060bt0aHTp0kDM/IiKiyqewqwkqPIHf3d0dr7/+Og4dOoQjR46gT58+MqRFRERkO0Iv3+YIZL2az8PDAwUFBXKGJCIiIiuTtRj4+++/sWLFCjlDEhERVT6FDRPIshwxERFRVeIo3fty4aJ/RERECseeASIiIlMK6xmQVAy88MILZbbfvHmzIrkQERHZBQ4TlMH0fgKmW5MmTRAWFmatXImIiKq8+Ph4+Pn5wdXVFQEBATh48GCZj7958yYmT56MBg0aQKvVomXLlti6daukc0rqGZg1axb8/PygVnOqARERVV226hlITk5GVFQUEhISEBAQgMWLFyMkJARnz56Fp6dnsccXFhaib9++8PT0xIYNG+Dj44NLly6hdu3aks4r6Vvd398f165dM+yPGDECOTk5kk5IRERk72y16NDChQsxceJEhIeHo23btkhISED16tWRmJhY4uMTExNx/fp1bNq0CT179oSfnx969+6Njh07SjqvpGJACOMbN2zduhX5+fmSTkhERGT3hEq2TafT4datW0abTqcrdsrCwkIcPnwYwcHBhmNqtRrBwcFIS0srMc3NmzcjMDAQkydPhpeXF9q3b4/58+ejqKhI0stlfz8REZEVxcXFFZtjFxcXV+xx165dQ1FREby8vIyOe3l5ITs7u8TY6enp2LBhA4qKirB161bMmjUL//nPfzBv3jxJOUqaM6BSqYrdrljq7YuJiIjsnZxzBqKjoxEVFWV0TKvVyhJbr9fD09MTn3/+OTQaDbp06YLLly/jo48+QmxsrMVxJBUDQgi8/PLLhhdRUFCA1157DTVq1DB63MaNG6WEJSIisitCL98fulqt1qIvfw8PD2g0mmJz8XJycuDt7V3icxo0aABnZ2doNBrDsTZt2iA7OxuFhYVwcXGxKEdJwwTjxo2Dp6enoZtjzJgxaNiwYbHuDyIiIpLGxcUFXbp0QWpqquGYXq9HamoqAgMDS3xOz549ceHCBej1/+vKOHfuHBo0aGBxIQBI7Bn44osvpDyciIjIIdnq0sKoqCiMGzcOXbt2Rbdu3bB48WLk5+cjPDwcABAWFgYfHx/DnIOIiAgsXboUb731Ft544w2cP38e8+fPx5tvvinpvFyOmIiIyIQQtpkPN2LECOTm5iImJgbZ2dno1KkTtm3bZphUmJWVZbTWj6+vL3788UdMmTIFHTp0gI+PD9566y1MnTpV0nlZDBAREdmRyMhIREZGlti2e/fuYscCAwNx4MCBCp2TxQAREZEJpd2bgMUAERGRCTmvJnAEXHSIiIhI4VTCdI1hIiIihcvq+oxssRofSjX/IBvjMAEREZEJpQ0T2FUx4OTiU6Hn3y+8bJdx7l1Lr1AcZ49mdpkP4zAO41gvjr19jtlLnEdjWZPSigHOGSAiIlI4u+oZICIisgdKm03HYoCIiMgEhwmIiIhIUdgzQEREZMJW9yawFRYDREREJpS2HDGHCYiIiBSOPQNEREQm9BwmICIiUjalzRngMAEREZHCsWeAiIjIhNLWGWAxQEREZIIrEBIRESmc0noGJM0Z+P7776HXK+ziSyIioipOUjEwZMgQ+Pr6YsaMGbhw4YK1ciIiIrIpvVDJtjkCScVARkYGJk2ahHXr1qFVq1bo3bs3Vq9ejbt371orPyIiokonhEq2zRFIKgZ8fX0RExODixcvYseOHfDz80NERAQaNGiA1157Db/++qu18iQiIiIrKfc6A0FBQVi5ciWuXr2Kjz76CCdOnED37t3RsWNHOfMjIiKqdELItzmCCl9NULNmTTzzzDO4dOkSzpw5g99//12OvIiIiGzGUcb65VLunoG7d+9i1apV6NOnD/z9/bFu3TpERUUhMzNTxvSIiIjI2iT3DBw4cACJiYlYv349CgsL8cILL2DHjh0ICgqyRn5ERESVzlEm/slFUjHQtm1bnD17Fp07d0ZcXBxeeukluLu7Wys3IiIim3CUsX65SCoGgoODsXbtWk4SJCIiqkIkFQNLliwp8fiePXuQn5+PwMBA1KlTx2wcnU4HnU5ndEyr1UpJhYiIyGo4gbAMH3zwAWbNmmXYF0Kgf//+CAoKwuDBg9GmTRucOnXKbJy4uDi4u7sbbXFxcdKzJyIisgIuOlSG5ORktG/f3rC/YcMG7N27F/v27cO1a9fQtWtXzJ4922yc6Oho5OXlGW3R0dHSsyciIrICpS1HLGmYICMjAx06dDDsb926FcOGDUPPnj0BADNnzsTw4cPNxtFqtRwWICIishOSegbu379v9CWelpaGHj16GPYbNmyIa9euyZcdERGRDQgZN0cgqRho3rw59u7dCwDIysrCuXPn0KtXL0P7n3/+iXr16smbIRERUSXjMEEZJk+ejMjISOzbtw8HDhxA9+7d0bZtW0P7zp070blzZ9mTJCIiIuuRVAxMnDgRGo0GKSkp6NWrF2JjY43ar1y5gvDwcFkTJCIiqmyOchWAXCQVA7du3cKwYcMwbNgwo2MPLViwQL7MiIiIbERv6wQqmaRioHbt2lCpzFdLRUVF5U6IiIiIKpekYmDXrl2G/xdCYODAgVi+fDl8fHxkT4yIiMhWBDhMUKrevXsb7Ws0GnTv3h3NmjWTNSkiIiJb0jvKNYEykXRpIREREVU9knoGiIiIlEDPYQJpLJlQSERE5Eg4Z6AML7zwgtF+QUEBXnvtNdSoUcPo+MaNGyueGRERkY3w0sIyuLu7G+2PGTNG1mSIiIio8kkqBr744gtr5UFERGQ3OExARESkcEobJuClhURERArHngEiIiITSusZYDFARERkQmlzBjhMQEREpHDsGSAiIjKhV1bHAIsBIiIiU0pbjpjDBERERAqnEkIo7EaNREREZdvk/ZJssYZkr5EtlrVwmICIiMgELy20IScXnwo9/37hZQDAvWvpFYrj7NFM1jh8XSV7+LoYp+w49vbzsrffQ3t7f/g+lx2nou8P8L/3yJr0CrsjL+cMEBERKZxd9QwQERHZA6VNpmMxQEREZEJpcwY4TEBERKRw7BkgIiIyobQVCNkzQEREZEIPlWybVPHx8fDz84OrqysCAgJw8OBBi563bt06qFQqDBkyRPI5WQwQERHZieTkZERFRSE2NhZHjhxBx44dERISgr/++qvM52VmZuL//u//8NRTT5XrvCwGiIiITAgZNykWLlyIiRMnIjw8HG3btkVCQgKqV6+OxMTEUp9TVFSE0aNHY/bs2WjWrJnEMz7AYoCIiMiEXiXfptPpcOvWLaNNp9MVO2dhYSEOHz6M4OBgwzG1Wo3g4GCkpaWVmuucOXPg6emJ8ePHl/v1shggIiKyori4OLi7uxttcXFxxR537do1FBUVwcvLy+i4l5cXsrOzS4z9888/Y8WKFVi2bFmFcuTVBERERCbkXGcgOjoaUVFRRse0Wm2F496+fRtjx47FsmXL4OHhUaFYLAaIiIhMyLkCoVartejL38PDAxqNBjk5OUbHc3Jy4O3tXezxFy9eRGZmJkJDQw3H9PoHZYyTkxPOnj2L5s2bW5QjhwmIiIhMyDlnwFIuLi7o0qULUlNT/5eHXo/U1FQEBgYWe3zr1q1x4sQJHDt2zLA9++yzCAoKwrFjx+Dr62vxudkzQEREZCeioqIwbtw4dO3aFd26dcPixYuRn5+P8PBwAEBYWBh8fHwQFxcHV1dXtG/f3uj5tWvXBoBix81hMUBERGTCVvcmGDFiBHJzcxETE4Ps7Gx06tQJ27ZtM0wqzMrKglotf6c+iwEiIiITtrxRUWRkJCIjI0ts2717d5nPTUpKKtc5OWeAiIhI4dgzQEREZEIo7EZFLAaIiIhM2HKYwBY4TEBERKRw7BkgIiIyobSeARYDREREJuRcgdARSBom2LlzJ9q2bYtbt24Va8vLy0O7du2wb98+2ZIjIiIi65NUDCxevBgTJ05ErVq1irW5u7tj0qRJWLhwoWzJERER2YItliO2JUnFwG+//Yb+/fuX2t6vXz8cPnzYbBxL7+1MRERkC3oZN0cgqRjIycmBs7Nzqe1OTk7Izc01G8fSezsTERHZAouBMvj4+ODkyZOlth8/fhwNGjQwGyc6Ohp5eXlGW3R0tJRUiIiISCaSioGBAwdi1qxZKCgoKNZ29+5dxMbGYvDgwWbjaLVa1KpVy2iz5F7PRERElUHIuDkCSZcWzpw5Exs3bkTLli0RGRmJVq1aAQDOnDmD+Ph4FBUVYcaMGVZJlIiIqLI4ysQ/uUgqBry8vLB//35EREQgOjoaQjyoeVQqFUJCQhAfH2+4zSIRERE5BsmLDjVp0gRbt27FjRs3cOHCBQgh4O/vjzp16lgjPyIiokrnKBP/5FLuFQjr1KmDJ554Qs5ciIiI7IKjjPXLhTcqIiIiUjjem4CIiMiEXmF9AywGiIiITChtzgCHCYiIiBSOPQNEREQmlDVIwGKAiIioGKUNE7AYICIiMqG0FQg5Z4CIiEjh2DNARERkgpcWEhERKZyySgEOExARESkeewaIiIhM8GoCIiIihVPanAEOExARESkcewaIiIhMKKtfgMUAERFRMZwzQEREpHCcM0BERESKohJCKKv8ISIiMmOK30jZYi3KXCdbLGvhMAEREZEJzhmwIScXnwo9/37hZQDAvWvpFYrj7NGsSudjb6+rqr4/fJ9LZq/vM+M4RpxHY5F87KoYICIisgdCYRMIWQwQERGZUNowAa8mICIiUjj2DBAREZlQ2joDLAaIiIhMKKsU4DABERGR4rFngIiIyASHCYiIiBROaVcTsBggIiIyobR1BjhngIiISOHYM0BERGSCwwREREQKx2ECIiIiUhT2DBAREZngMAEREZHC6QWHCYiIiEhB2DNARERkQln9AiwGiIiIilHacsQcJiAiIlI49gwQERGZUNo6AywGiIiITPDSQiIiIoXjnIEypKWl4fvvvzc6tmrVKjRt2hSenp549dVXodPpZE2QiIiIrEtSMTBnzhycOnXKsH/ixAmMHz8ewcHBmDZtGlJSUhAXFyd7kkRERJVJyPifI5A0THDs2DHMnTvXsL9u3ToEBARg2bJlAABfX1/ExsbivffeKzOOTqcr1oOg1WqlpEJERGQ1SpszIKln4MaNG/Dy8jLs79mzBwMGDDDsP/HEE/jjjz/MxomLi4O7u7vRxh4FIiIi25BUDHh5eSEjIwMAUFhYiCNHjqB79+6G9tu3b8PZ2dlsnOjoaOTl5Rlt0dHRElMnIiKyDiGEbJtU8fHx8PPzg6urKwICAnDw4MFSH7ts2TI89dRTqFOnDurUqYPg4OAyH18aScXAwIEDMW3aNOzbtw/R0dGoXr06nnrqKUP78ePH0bx5c7NxtFotatWqZbRxmICIiOyFHkK2TYrk5GRERUUhNjYWR44cQceOHRESEoK//vqrxMfv3r0bo0aNwq5du5CWlgZfX1/069cPly9flnReScXA3Llz4eTkhN69e2PZsmVYtmwZXFxcDO2JiYno16+fpASIiIjogYULF2LixIkIDw9H27ZtkZCQgOrVqyMxMbHEx3/11Vd4/fXX0alTJ7Ru3RrLly+HXq9HamqqpPNKmkDo4eGBvXv3Ii8vD25ubtBoNEbtX3/9Ndzc3CQlQEREZG/knEBY2qR50x7xwsJCHD582GjYXK1WIzg4GGlpaRad686dO7h37x7q1q0rKcdy3ZvA3d29WCEAAHXr1jXqKSAiInJEcl5aaOmk+WvXrqGoqMhooj7wYL5edna2RXlPnToVDRs2RHBwsKTXyxUIiYiIrCg6OhpRUVFGx6wxT27BggVYt24ddu/eDVdXV0nPZTFARERkQs7liEsaEiiJh4cHNBoNcnJyjI7n5OTA29u7zOf++9//xoIFC7Bjxw506NBBco68hTEREZEJW1xa6OLigi5duhhN/ns4GTAwMLDU53344YeYO3cutm3bhq5du5br9bJngIiIyIStViCMiorCuHHj0LVrV3Tr1g2LFy9Gfn4+wsPDAQBhYWHw8fExzDn44IMPEBMTgzVr1sDPz88wt8DNzU3ShH4WA0RERHZixIgRyM3NRUxMDLKzs9GpUyds27bNMKkwKysLavX/OvU/++wzFBYWYtiwYUZxLLk1wKNYDBAREZmw5Q2GIiMjERkZWWLb7t27jfYzMzNlOSeLASIiIhNyTiB0BJxASEREpHDsGSAiIjJRnhsMOTIWA0RERCY4TEBERESKwp4BIiIiE7a8msAWWAwQERGZ0CtszgCHCYiIiBSOPQNEREQmlNUvwGKAiIioGKVdTcBigIiIyITSigHOGSAiIlI49gwQERGZUNoKhCqhtFdMRERkRreGvWWLdfDKHtliWQuHCYiIiBTOroYJnFx8KvT8+4WX7TLOvWvpFYrj7NGsSudjb6/L3uLY2/vMOJUTx95+D+0lzqOxrIkrEBIRESmc0kbQOUxARESkcOwZICIiMqG0dQZYDBAREZngMAEREREpCnsGiIiITHCYgIiISOF4aSEREZHC6TlngIiIiJSEPQNEREQmlDZMUK6egSNHjuDEiROG/e+++w5DhgzB9OnTUVhYKFtyREREtqAXQrbNEZSrGJg0aRLOnTsHAEhPT8fIkSNRvXp1fP3113j33XdlTZCIiIisq1zFwLlz59CpUycAwNdff41evXphzZo1SEpKwjfffCNnfkRERJVOyPifIyjXnAEhBPR6PQBgx44dGDx4MADA19cX165dky87IiIiG3CU7n25lKtnoGvXrpg3bx5Wr16NPXv2YNCgQQCAjIwMeHl5yZogERERWVe5ioHFixfjyJEjiIyMxIwZM9CiRQsAwIYNG9CjRw9ZEyQiIqpsHCawQIcOHYyuJnjoo48+gkajqXBSREREtqS0YQJZ1xlwdXWVMxwRERFVgnIVA0VFRVi0aBHWr1+PrKysYmsLXL9+XZbkiIiIbMFRuvflUq45A7Nnz8bChQsxYsQI5OXlISoqCi+88ALUajXee+89mVMkIiKqXELoZdscQbmKga+++grLli3Dv/71Lzg5OWHUqFFYvnw5YmJicODAAblzJCIiqlR6CNk2R1CuYiA7OxuPPfYYAMDNzQ15eXkAgMGDB2PLli3yZUdERERWV65ioFGjRrh69SoAoHnz5ti+fTsA4Ndff4VWq5UvOyIiIhsQQsi2OYJyTSB8/vnnkZqaioCAALzxxhsYM2YMVqxYgaysLEyZMsXs83U6HXQ6ndExFhFERGQvHKV7Xy7lKgYWLFhg+P8RI0agcePGSEtLg7+/P0JDQ80+Py4uDrNnzzY6FhsbW55UiIiIqIJkWWcgMDAQgYGBFj8+OjoaUVFRRse0Wi3mzV8mRzpEREQV4ijd+3KxuBjYvHmzxUGfffbZMtu1Wi2HBYiIyG5xBcJSDBkyxKLHqVQqFBUVlTcfIiIiqmQWFwMPb1lMRERU1SltBUJZ701ARERUFXDOQCmWLFlicdA333yzXMkQERFR5bO4GFi0aJFFj1OpVCwGiIjIoXGdgVJkZGQUO5abmwuVSgUPDw9ZkyIiIrIlpQ0TSF6O+ObNm5g8eTI8PDzg7e0NLy8veHh4IDIy0nCPAiIiIkemF0K2zRFImkB4/fp1BAYG4vLlyxg9ejTatGkDAPj999+RlJSE1NRU7N+/H3Xq1LFKskRERCQ/ScXAnDlz4OLigosXL8LLy6tYW79+/TBnzhyL5xcQERHZIw4TlGHTpk3497//XawQAABvb298+OGH+Pbbb2VLjoiIyBb0ELJtjkBSMXD16lW0a9eu1Pb27dsjOzu7wkkRERFR5ZFUDHh4eCAzM7PU9oyMDNStW7eiOREREdmUEEK2zRFIKgZCQkIwY8YMFBYWFmvT6XSYNWsW+vfvL1tyREREtsCrCcowZ84cdO3aFf7+/pg8eTJat24NIQROnz6NTz/9FDqdDqtXr7ZWrkRERGQFkoqBRo0aIS0tDa+//jqio6MN3R8qlQp9+/bF0qVL4evra5VEiYiIKgtvVGRG06ZN8cMPP+DGjRs4f/48AKBFixacK0BERFWGo3Tvy6Xcdy2sU6cOunXrJmcuREREZAO8hTEREZEJR7kKQC4sBoiIiExwzgAREZHCKa1nQPJdC4mIiMh64uPj4efnB1dXVwQEBODgwYNlPv7rr79G69at4erqisceewxbt26VfE4WA0RERCZstQJhcnIyoqKiEBsbiyNHjqBjx44ICQnBX3/9VeLj9+/fj1GjRmH8+PE4evQohgwZgiFDhuDkyZOSzstigIiIyISQcZNi4cKFmDhxIsLDw9G2bVskJCSgevXqSExMLPHxH3/8Mfr374933nkHbdq0wdy5c/H4449j6dKlks6rEkobGCEiIjLDycVHtlj5t9Oh0+mMjmm1Wmi1WqNjhYWFqF69OjZs2IAhQ4YYjo8bNw43b97Ed999Vyx248aNERUVhbfffttwLDY2Fps2bcJvv/1mcY4O0zOg0+nw3nvvFXtDGYdxGIdxGKdqxJE7VkXcL7ws2xYXFwd3d3ejLS4urtg5r127hqKiInh5eRkd9/LyKvWOwNnZ2ZIeXyrhIPLy8gQAkZeXxziMwziMwzhVMI7csexFQUGByMvLM9oKCgqKPe7y5csCgNi/f7/R8XfeeUd069atxNjOzs5izZo1Rsfi4+OFp6enpBx5aSEREZEVlTQkUBIPDw9oNBrk5OQYHc/JyYG3t3eJz/H29pb0+NI4zDABERFRVebi4oIuXbogNTXVcEyv1yM1NRWBgYElPicwMNDo8QDw008/lfr40rBngIiIyE5ERUVh3Lhx6Nq1K7p164bFixcjPz8f4eHhAICwsDD4+PgY5hy89dZb6N27N/7zn/9g0KBBWLduHQ4dOoTPP/9c0nkdphjQarWIjY21qKuFcRiHcRiHcRwvjtyxHNGIESOQm5uLmJgYZGdno1OnTti2bZthkmBWVhbU6v916vfo0QNr1qzBzJkzMX36dPj7+2PTpk1o3769pPPy0kIiIiKF45wBIiIihWMxQEREpHAsBoiIiBSOxQAREZHCsRggIiJSOIe5tNAeZGdn47///a9hzWdvb28EBARIXunJGjIyMnDhwgU0aNBA8iUlRKQMv/32Gw4fPow+ffqgWbNmOHXqFOLj46HX6/H8888jJCTE4lg7d+7Ezz//jKtXr0KtVqNZs2Z49tln4e/vb8VXQFYjafHiShIfHy+eeeYZMXz4cLFjxw6jttzcXNG0aVOzMf773/+K+/fvG/ZTUlJEr169RMOGDUWXLl3EypUrLc7nn3/+EaNHjxYajUY4OTkJT09P4enpKZycnIRGoxFjxowR+fn5ZuMUFhaKd955RzRv3lw88cQTYsWKFUbt2dnZQq1Wm40TEREhbt++LYQQ4s6dO2Lo0KFCrVYLlUol1Gq1CAoKMrRXxPXr1yW9T1evXhWbNm0SCQkJIiEhQWzatElcvXq1wnk89M8//4g9e/ZY9Fi9Xi/S09PFvXv3hBBC6HQ6sW7dOrFy5UqRm5tboTyCgoJEZmamxY/fsGGDRb8fljp27JhYsWKFuHjxohBCiJMnT4qIiAgxadIksW3bNkmxioqKSj1+6dKlcuco9T0qSXp6uti+fbs4ceKEpOfJ+f7YA7k+N7755huh0WhEvXr1hJubm/jpp59E7dq1RXBwsAgJCREajUZ89dVXZuPk5OSIbt26CbVaLZycnIRarRZdunQR3t7eQqPRiHfeeafcr1UIIc6dOyd27Nghzp8/X6E4JI3dFQMff/yxqF69upg8ebIYM2aMcHFxEfPnzze0W/qLr1arRU5OjhBCiM2bNwu1Wi3CwsJEfHy8mDBhgnBychIbN260KKfx48cLf39/sW3bNqMC4/79++LHH38ULVu2FBMmTDAbJzY2Vnh5eYmPPvpIzJgxQ7i7u4tXX33V6LWpVCpJry06Olo0atRI7Ny5U+Tn54uff/5ZNG/eXEybNs2i11aWY8eOWfRey1UsyZXPmTNnRJMmTYRarRYtWrQQ6enpokuXLqJGjRqievXqwsPDQ5w7d85snO+++67ETaPRiKVLlxr2zVGpVKJWrVpi4sSJ4sCBAxa91tLI9YGel5cnhg8fLlxdXYWnp6eYNWuW0e+2pf/O5HqP5Cpw5Xp/TP3zzz8iMTFRTJ8+XXzyySfi2rVrFj/3k08+EWPHjhVr164VQgixatUq0aZNG9GqVSsRHR1tKFhLI9fnxuOPPy7mzZsnhBBi7dq1onbt2mLOnDmG9n//+9+iU6dOZuOMGDFCDBkyxHCzncjISBEWFiaEECI1NVXUq1dPLF682GwcIYSYP3++4Q++69evi2eeeUaoVCrDz71///7ixo0bFsWiirG7YqBt27ZG/1h/+eUXUb9+fTFr1iwhhOUfUiqVyvCF+eSTTxb7cnz//fdF9+7dLcqpdu3a4pdffim1/eeffxa1a9c2G6dFixYiJSXFsH/+/HnRokUL8fLLLwu9Xl+u19a+fftid6z67rvvRMuWLc3GMb2Llum2b98+i/KRq1gyx9Ji4LnnnhPPPvusOH78uHj77bdFmzZtxHPPPScKCwtFQUGBCA0NFWPGjDEb5+EH0sMPp5I2S39ec+bMEZ07dxYqlUq0a9dOLFq0SNIXykNyfaC/+eabomXLluLrr78Wy5YtE02aNBGDBg0SOp1OCGH5F4xc75FcBa5c70+bNm3E33//LYQQIisrS/j5+Ql3d3fxxBNPiLp16wpPT0+Rnp5uNs7cuXNFzZo1xdChQ4W3t7dYsGCBqFevnpg3b56YP3++qF+/voiJiSkzhlyfGzVq1BAZGRlCiAc9Z87OzuL48eOG9osXLwo3NzezcWrVqiVOnjxp2P/nn3+Es7Oz4S6Dq1evFq1atTIbRwghGjVqJI4cOSKEEGLChAmic+fO4siRI+Lu3bvi2LFjonv37mL8+PEWxaKKsbtioFq1aoZf2IdOnDghvLy8xLRp08r1henp6SkOHTpk1H7mzBmLvsCFePDL/+uvv5bafvDgQVGrVi2zcUp6bX/++ado2bKlGD16tLh8+bLFr+2vv/4SQgjh4eFh9A9TCCEyMzNFtWrVLIqjVqtL3Sz9IJerWKpTp06ZW61atSzKp379+uLo0aNCiAcfVCqVSuzbt8/Q/ssvv4jGjRubjdO/f38xaNAgw+/RQ05OTuLUqVNmn//Qo7+Lhw4dEhEREaJ27dpCq9WK4cOHi+3bt1scS64P9MaNG4tdu3YZ9nNzc0W3bt1Ev379REFBgcX/zqzxHlWkwJXr/Xk0n9GjR4sePXqImzdvCiGEuH37tggODhajRo0yG6d58+bim2++EUI8KGY1Go348ssvDe0bN24ULVq0KDOGXJ8b3t7ehs/B69evC5VKZfQ7cPDgQeHt7W02Tv369Y1+tnfu3BFqtdpQPF28eFFotVqzcYQQQqvVGoaT/Pz8ig0DHjp0SDRo0MCiWFQxdjeB0MPDA3/88Qf8/PwMx9q3b4+dO3fi6aefxpUrVyyO9fvvvyM7OxvVqlWDXq8v1n7//n2L4gwePBivvvoqVqxYgc6dOxu1HT16FBEREQgNDTUbx9vbGxcvXjR6bT4+Pti1axeCgoLw8ssvW5QPAMyaNQvVq1eHWq3GlStX0K5dO0Pb33//jRo1apiNUbNmTcyYMQMBAQEltp8/fx6TJk0yG0ev18PFxaXUdhcXlxLff1M6nQ4RERF47LHHSmy/dOkSZs+ebTbOP//8g7p16wIAatSogRo1aqBBgwaGdl9f32K3/CzJDz/8gEWLFqFr16749NNPMXjwYLPPMadLly7o0qULFi5ciK+//hqJiYno378/GjdujIyMDLPPr1mzJv7++2/4+fnh5s2buH//Pv7++29D+99//w03NzezcXJzc9GkSRPDvoeHB3bs2IGQkBAMHDgQy5cvt+j1yPkeqVQqAA8m6nbo0MGorWPHjvjjjz/MxpDr/XlUWloaEhIS4O7uDgBwc3PD7NmzMXLkSLPPvXLlCrp27Wp4DWq1Gp06dTK0P/7442Y/0+T63AgODsbkyZPxxhtvIDk5Gf369UN0dDS++OILqFQqvPPOO3jyySfNxnnyyScRExODlStXwsXFBdOnT0ezZs0M/+Zyc3NRp04di3Jq0qQJTp48iSZNmkClUsHJyfgrSaPRID8/36JYVEG2rkZMjRo1Srz99tsltp08eVLUr1/f4r+eH+2+XLRokVH72rVrRdu2bS3K6fr166J///5CpVKJunXritatW4vWrVuLunXrCrVaLQYMGGDRuNb48ePFK6+8UmLbn3/+KVq0aGHRa+vdu7fo06ePYVu2bJlR+9y5c0Xv3r3NxunTp4/44IMPSm0/duyYRV3FL730kqF7z9SRI0dEly5dxOjRo83G6dGjR5ljjZYOEzRv3tyoJ+DTTz8Vt27dMuwfPnzYor+AHjp69Kho27atePXVV0V+fr7kv3of7QIvyfnz58X06dMtijVmzBgREBAgvvzySxEaGipCQkJE9+7dxenTp8WZM2dE7969xbBhw8zGadWqldiyZUux47dv3xaBgYGiY8eOFr3XD1X0PVKpVGLSpEliypQpwtPTs1hvyeHDh4WHh4fZOHK9P4/2vjVs2LDYJMbMzEzh6upqNk7Tpk3FDz/8IIR4MDFOrVaL9evXG9q3bNki/Pz8yowh1+dGdna26Nu3r3BzcxMhISHi5s2bIjIy0vAZ6e/vLy5cuGA2zsWLF0WzZs2Ek5OTcHZ2Fu7u7kY/ry+++MLiOUsfffSRaNOmjTh//rz4z3/+IwIDAw05pKeniz59+lj086KKs7ti4LfffhOJiYmltp88eVIMHTrUbJzMzEyjzXR8duXKlZJmygshxOnTp0ViYqKYP3++mD9/vkhMTBSnT5+2+PmZmZllzma+fPmySEpKkpRTSdLT08Uff/xh9nGff/65+Pjjj0ttz87OFu+9957ZOHIVS++//36Z58vKyhIvv/yy2TiTJk0qViA9Ki4uTgwcONBsnEfl5+eLV199Vfj7+wuNRlPuLvCKKusDXa1WW/yBHhkZWeqH7K1bt0RAQICkYkCIB93FkyZNKtd7JFeBK9f7o1KpxGOPPSY6d+4s3NzcxIYNG4za9+zZI3x8fMzGmTlzpqhfv76YMGGCaNq0qZg2bZpo3Lix+Oyzz0RCQoLw9fUVU6ZMKTOGtT83Ll68KE6cOGF2IuOj8vPzxfbt20VKSkqFr8554403hLOzs2jdurVwdXUVarVauLi4GK5SkPOKJCqdw9y18Pbt21i7di2WL1+Ow4cPo6ioqMIxr1+/bujaciQ7d+5EZGQkDhw4gFq1ahm15eXloUePHkhISMBTTz1VqXmdOXMGaWlpRuswBAYGonXr1pWahzkZGRlwdXU1GjqwVEpKCnbu3Ino6Gh4enpa9JxLly7B19fX6LajcktPT8edO3fQunXrYl2tJblx40ax4aVH3b59G0eOHEHv3r0l57J582bs2rVL0ntkTnp6OrRaLXx8fMr9fCnvj+lwVEBAAPr372/Yf+edd/Dnn39i7dq1ZcbR6/VYsGAB0tLS0KNHD0ybNg3Jycl49913cefOHYSGhmLp0qUWDetV1N27d5GammoYyomOjoZOpzO0Ozk5Yc6cOXB1da2UOI86ffo0vv/+e6Snp0Ov16NBgwbo2bMngoODDcNHZGW2rkbM2bNnjwgLCxM1atQQ/v7+YurUqeLgwYMVivnjjz+KF1980aJuPktIvR6/onFCQ0PFwoULS23/+OOPxZAhQyw6Z1FRkVixYoUYNGiQaNeunWjfvr0IDQ0VK1euFHq93uLc7Y1cr8ve4sgVa8CAAYYJcUI86C15tPfm2rVrok2bNoqNc/HixVLXYJBCrjhlsfRz47PPPhODBw827Lu5uYmAgABDb4y3t3eZnytyxxHC/BVNDzeyPrssBq5evSri4uJEixYthKenp4iMjJQ8BmkqMzNTxMTEiCZNmohatWqJESNGGI3dVYSlY9lyxWncuLH4/fffS20/ffq08PX1NRtHr9eLgQMHCpVKJTp16iRGjhwpRowYITp06CBUKpV47rnnpKRfqsouluR6XXq9XgwaNMhu4sj52kznMdSsWdOwSI8Q5VvPoyJxTIdSbJ2PaZwXX3xRZGdnm32eteKUxdLPjSeffFJs3rzZsO/m5mb03qxevdqiy63liiOEfFc0UcXZXTEwePBgUatWLTFq1Cjx/fffG65bL08xoNPpxNq1a8UzzzwjXF1dxeDBg4VGozG61MgScl2PL1ccrVZb5upc58+ft6jXIzExUdSsWVPs3LmzWFtqaqqoWbOmLF/ilV0syfW67C2OnLFMv3xNP9DL+yWulDiWkiOOXJ8b3t7eRpcoenh4GO2fPXvWokuk5YojhBC7d+82bLt27RLVqlUTX331ldHx3bt3WxSLKsbuigGNRiOmTJlSbIU4qcVAZGSkqFevnujevbtYunSpYQJheYoKuapXueI0a9ZMfPvtt6W2f/PNNxYt2dy3b18RFxdXavv7778v+vXrZzaOvRVLcr0ue4sjZyxH+/K1tziWkiOOXJ8brq6u4syZM6W2nz592qL1AeSKU5Lyvs9UcXa3zsDPP/+MFStWoEuXLmjTpg3Gjh1r0fW8pj777DNMnToV06ZNQ82aNSuUk1zX48sVZ+DAgZg1axb69+9fbJLO3bt3ERsba9H13sePH8eHH35YavuAAQOwZMkSs3Fq165d5iQfIYRFk4DkiiPX67K3OHLGUqlUxd7L8kzUYhzrx5Hrc6NRo0Y4efIkWrVqVWL78ePH0ahRo0qLQ/bF7oqB7t27o3v37li8eDGSk5ORmJiIqKgo6PV6/PTTT/D19bXoy3316tVITExEgwYNMGjQIIwdOxYDBgwoV06PP/44AJQ6s7p27doQFlyUIVecmTNnYuPGjWjZsiUiIyMN/yjPnDmD+Ph4FBUVYcaMGWbjXL9+HV5eXqW2e3l54caNG2bj2FuxJNfrsrc4csYSQuDll1+GVqsFABQUFOC1114zzGp/dHY44xSP89DGjRutHkeuz42BAwciJiYGgwYNKvGPiNmzZ2PQoEGVFofsi90VAw/VqFEDr7zyCl555RWcPXsWK1aswIIFCzBt2jT07dsXmzdvLvP5o0aNwqhRo5CRkYGkpCRMnjwZd+7cgV6vx++//462bdtanMtLL72EO3fulNru7e2N2NjYSovj5eWF/fv3IyIiAtHR0YYPApVKhZCQEMTHx5f5hfFQUVFRmZdZaTQai1ZptLdiSa7XZW9x5Iw1btw4o/0xY8YUe0xYWBjjlBHHEnLEketzY/r06Vi/fj1atWqFyMhItGzZEgBw9uxZLF26FPfv38f06dMrLU5peCmhbTjMOgPAgw/ClJQUJCYmmi0GTAkhsH37dqxYsQKbN2+Gh4cHXnjhBYu7Z+3VjRs3cOHCBQgh4O/vb/EyoACgVqsxYMAAw18tpnQ6HbZt22Z2TYdly5bh7t27ePPNN0tsz8nJQUJCgtkPLLniyPW67C2O3LFIeTIyMhAREYGffvrJ6I+Ivn374tNPP0WzZs0qNc4LL7xgtJ+SkoKnn35acg8MVZxDFQNyuX79OlavXo1//etfFv0VJdciP/a2WFB4eLhFj/viiy+snIm85Hpd9hZH7ljkGAoKCrBjx45SF/nRaDSYO3eupEV+rl+/jgsXLgAAWrRoUe7F1yoah7/P9kNxxUB5VjJ89tlnERQUhClTppTYvmTJEuzatQvffvttpcSxR3q9HklJSdi4cSMyMzOhUqnQtGlTDBs2DGPHjrW460+uOERVRUJCArZs2YKUlBQAD+bWtGvXDtWqVQPwYK7Qu+++W+rnCpFFKu/CBduqyEqGci3yI1cceyPnIj+VsQgSkSORc5EfotLY7QRCOWRnZyMpKQkrVqzArVu38OKLL0Kn02HTpk2SJhDm5OTA2dm51HYnJyfk5uZWWhx7k5SUhH379iE1NRVBQUFGbTt37sSQIUOwatUqsxO35IpDVJVcuHDB6Lberq6uRve56NatGyZPnmyL1KgKsd6dU2wsNDQUrVq1wvHjx7F48WJcuXIFn3zySbli+fj44OTJk6W2Hz9+3KKb3sgVx96sXbsW06dPL/YFDgBPP/00pk2bhq+++qrS4hBVJTdv3jSaI5Cbmws/Pz/Dvl6vt/iSSaLSVNli4IcffsD48eMN17xqNJpyx3q4yE9BQUGxNimL/MgVx94cP37c6I5upgYMGIDffvut0uIQVSUPF/kpDRf5ITlU2QmEBw4cwIoVK5CcnGy0kmGDBg3w22+/SR4mePzxx6HRaEpd5OfIkSNmr+2XK469cXFxwaVLl0rt1bhy5QqaNm1q9q8XueIQVSVvvfUWduzYgcOHD5e4yE/Xrl0RHByMjz/+2EYZUlVQZYuBh/Lz8w0rGR48eBBFRUVYuHAhXnnlFUnLFF+6dAkRERH48ccfS1zkp2nTppUax55oNBpkZ2ejfv36Jbbn5OSgYcOGZq/ckCsOUVWSk5ODTp06wcXFpdRFfo4ePepwf0SQfanyxcCjHq5kuHr1aty8edOilQxNVWSRH2vEsQf2uDgPUVUi1yI/RKVRVDHwUEVWMqTi7HFxHqKqSK7FgohMKbIYICIiov+pslcTEBERkWVYDBARESkciwEiIiKFYzFARESkcCwGiIiIFI7FABERkcKxGCAiIlI4FgNEREQK9//xsZQ/yYqPXgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# define inputs\n",
    "##########################################################################################\n",
    "IN_DIR = 'merged'\n",
    "OUT_DIR = 'mask'\n",
    "META_DIR = 'metadata'\n",
    "# define input directory\n",
    "#DATA_DIR = '101222_D10_Coverslip1_Processed'\n",
    "##########################################################################################\n",
    "# os.chdir(f'{DATA_DIR}')\n",
    "# print(os.getcwd())\n",
    "\n",
    "#SOURCE = f'gs://fc-secure-9289bfef-e5cb-493a-83d5-e604cd429e39/Brian/{DATA_DIR}/'\n",
    "\n",
    "# load metadata - load full codebook as well\n",
    "full_codebook = pd.read_csv(f'{META_DIR}/full_codebook.csv',sep=',', index_col=0) # this is \"legal\" codebook\n",
    "Procode_gRNA = pd.read_csv(f'{META_DIR}/PROCODE_gRNA.csv',sep=',')\n",
    "legal_codes = sorted(list(set(Procode_gRNA['ProCode ID'].to_list())))\n",
    "codebook = full_codebook[legal_codes]\n",
    "#AllProcodes = pd.read_csv('AllProcodes.csv', sep='.')\n",
    "sns.heatmap(codebook, linewidths = 0.3)\n",
    "markers = pd.read_csv(f'{META_DIR}/markers.csv')\n",
    "#markers.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0,\n",
       " 5,\n",
       " 6,\n",
       " 7,\n",
       " 9,\n",
       " 10,\n",
       " 11,\n",
       " 13,\n",
       " 14,\n",
       " 15,\n",
       " 17,\n",
       " 18,\n",
       " 19,\n",
       " 21,\n",
       " 22,\n",
       " 23,\n",
       " 25,\n",
       " 27,\n",
       " 29,\n",
       " 30,\n",
       " 31,\n",
       " 33,\n",
       " 35,\n",
       " 36,\n",
       " 37]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# get indices of final channels --> only make mask when happening in these indices of merged file\n",
    "final_inds = []\n",
    "marker_names = ['DNA_0','NWS','VSVG','FLAG','HSV','C','S','Ollas','GFAP','NeuN',\n",
    "               'pRPS6','RANGAP1','NFKB','TOM20','LAMP1','4HNE','TDP43','G3BP1','GM130','Calnexin','Golgin97',\n",
    "               'SYTO','ER','AGP','Catalase']\n",
    "\n",
    "for marker in marker_names:\n",
    "    ind_to_add = markers[markers.marker_name == marker]['Split_Num'].values[0]\n",
    "    final_inds.append(ind_to_add)\n",
    "final_inds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['001', '002', '004', '005', '006', '007', '008', '009', '010', '012', '014', '015', '016', '017', '019', '020', '021', '022', '023', '024', '026', '027', '029', '030', '031', '032', '033', '034', '035', '036', '038', '039', '040', '041', '042', '043', '044', '045', '046', '047', '048', '049', '051', '053', '054', '055', '056', '057', '058', '059', '060', '061', '063', '064', '065', '066', '067', '068', '069', '070', '071', '072', '073', '074', '075', '076', '078', '079', '080', '081', '083', '084', '085', '086', '087', '088', '089', '090', '091', '092', '093', '094', '095', '096', '097', '100', '101', '102', '103', '104', '105', '106', '107', '108', '109', '110', '111', '112', '113', '114', '115', '116', '117', '118', '119', '120', '121', '122', '123', '124', '125', '126', '127', '128', '129', '130', '131', '132', '133', '134', '135', '136', '137', '138', '139', '140', '141', '142', '143', '144', '145', '146', '147', '148', '149', '150', '151', '152', '153', '154', '155', '156', '158', '159', '160', '161', '162', '163', '164', '165', '166', '167', '168', '169', '170', '171', '172', '173', '174', '175', '176', '177', '178', '179', '180', '181', '182', '183', '184', '185', '186', '188', '189', '190', '191', '192', '193', '194', '195', '196', '197', '199', '201', '202', '207', '208', '209', '210', '211', '212', '217', '218', '219', '220', '221']\n",
      "195\n"
     ]
    }
   ],
   "source": [
    "_allFOVs = sorted(glob.glob('tmat_Cyc_2/*'))\n",
    "allFOVs = [x.split('F')[-1][:3] for x in _allFOVs]\n",
    "print(allFOVs)\n",
    "print(len(allFOVs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['001', '002', '004', '005', '006', '007', '008', '009', '010', '012']\n",
      "['014', '015', '016', '017', '019', '020', '021', '022', '023', '024']\n",
      "['026', '027', '029', '030', '031', '032', '033', '034', '035', '036']\n",
      "['038', '039', '040', '041', '042', '043', '044', '045', '046', '047']\n",
      "['048', '049', '051', '053', '054', '055', '056', '057', '058', '059']\n",
      "['060', '061', '063', '064', '065', '066', '067', '068', '069', '070']\n",
      "['071', '072', '073', '074', '075', '076', '078', '079', '080', '081']\n",
      "['083', '084', '085', '086', '087', '088', '089', '090', '091', '092']\n",
      "['093', '094', '095', '096', '097', '100', '101', '102', '103', '104']\n",
      "['105', '106', '107', '108', '109', '110', '111', '112', '113', '114']\n",
      "['115', '116', '117', '118', '119', '120', '121', '122', '123', '124']\n",
      "['125', '126', '127', '128', '129', '130', '131', '132', '133', '134']\n",
      "['135', '136', '137', '138', '139', '140', '141', '142', '143', '144']\n",
      "['145', '146', '147', '148', '149', '150', '151', '152', '153', '154']\n",
      "['155', '156', '158', '159', '160', '161', '162', '163', '164', '165']\n",
      "['166', '167', '168', '169', '170', '171', '172', '173', '174', '175']\n",
      "['176', '177', '178', '179', '180', '181', '182', '183', '184', '185']\n",
      "['186', '188', '189', '190', '191', '192', '193', '194', '195', '196']\n",
      "['197', '199', '201', '202', '207', '208', '209', '210', '211', '212']\n",
      "['217', '218', '219', '220', '221']\n"
     ]
    }
   ],
   "source": [
    "for ii in range(0,len(allFOVs),10): # For each CHUNK (10 FOVs per loop)\n",
    "    print(allFOVs[ii:ii+10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "      <th>...</th>\n",
       "      <th>28</th>\n",
       "      <th>29</th>\n",
       "      <th>30</th>\n",
       "      <th>31</th>\n",
       "      <th>32</th>\n",
       "      <th>33</th>\n",
       "      <th>34</th>\n",
       "      <th>35</th>\n",
       "      <th>36</th>\n",
       "      <th>37</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>001</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>002</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>004</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>005</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>006</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>217</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>218</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>219</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>220</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>221</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>195 rows × 38 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      0    1    2    3    4    5    6    7    8    9   ...   28   29   30  \\\n",
       "001  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  ...  NaN  NaN  NaN   \n",
       "002  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  ...  NaN  NaN  NaN   \n",
       "004  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  ...  NaN  NaN  NaN   \n",
       "005  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  ...  NaN  NaN  NaN   \n",
       "006  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  ...  NaN  NaN  NaN   \n",
       "..   ...  ...  ...  ...  ...  ...  ...  ...  ...  ...  ...  ...  ...  ...   \n",
       "217  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  ...  NaN  NaN  NaN   \n",
       "218  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  ...  NaN  NaN  NaN   \n",
       "219  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  ...  NaN  NaN  NaN   \n",
       "220  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  ...  NaN  NaN  NaN   \n",
       "221  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  ...  NaN  NaN  NaN   \n",
       "\n",
       "      31   32   33   34   35   36   37  \n",
       "001  NaN  NaN  NaN  NaN  NaN  NaN  NaN  \n",
       "002  NaN  NaN  NaN  NaN  NaN  NaN  NaN  \n",
       "004  NaN  NaN  NaN  NaN  NaN  NaN  NaN  \n",
       "005  NaN  NaN  NaN  NaN  NaN  NaN  NaN  \n",
       "006  NaN  NaN  NaN  NaN  NaN  NaN  NaN  \n",
       "..   ...  ...  ...  ...  ...  ...  ...  \n",
       "217  NaN  NaN  NaN  NaN  NaN  NaN  NaN  \n",
       "218  NaN  NaN  NaN  NaN  NaN  NaN  NaN  \n",
       "219  NaN  NaN  NaN  NaN  NaN  NaN  NaN  \n",
       "220  NaN  NaN  NaN  NaN  NaN  NaN  NaN  \n",
       "221  NaN  NaN  NaN  NaN  NaN  NaN  NaN  \n",
       "\n",
       "[195 rows x 38 columns]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# record in dataframe\n",
    "count0_df = pd.DataFrame(index=allFOVs, columns = list(range(38)))\n",
    "count20_df = pd.DataFrame(index=allFOVs, columns = list(range(38)))\n",
    "count65K_df = pd.DataFrame(index=allFOVs, columns = list(range(38)))\n",
    "count65K_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FOV 001\n",
      "mask shape (2008, 2000)\n",
      "F001, no mask needed\n",
      "FOV 002\n",
      "mask shape (2008, 2000)\n",
      "F002, no mask needed\n",
      "FOV 004\n",
      "mask shape (2008, 2000)\n",
      "F004, no mask needed\n",
      "FOV 005\n",
      "mask shape (2008, 2000)\n",
      "F005, no mask needed\n",
      "FOV 006\n",
      "mask shape (2008, 2000)\n",
      "F006, no mask needed\n",
      "FOV 007\n",
      "mask shape (2008, 2000)\n",
      "F007, no mask needed\n",
      "FOV 008\n",
      "mask shape (2008, 2000)\n",
      "FOV 009\n",
      "mask shape (2008, 2000)\n",
      "FOV 010\n",
      "mask shape (2008, 2000)\n",
      "F010, no mask needed\n",
      "FOV 012\n",
      "mask shape (2008, 2000)\n",
      "F012, no mask needed\n",
      "FOV 014\n",
      "mask shape (2008, 2000)\n",
      "FOV 015\n",
      "mask shape (2008, 2000)\n",
      "F015, no mask needed\n",
      "FOV 016\n",
      "mask shape (2008, 2000)\n",
      "F016, no mask needed\n",
      "FOV 017\n",
      "mask shape (2008, 2000)\n",
      "F017, no mask needed\n",
      "FOV 019\n",
      "mask shape (2008, 2000)\n",
      "F019, no mask needed\n",
      "FOV 020\n",
      "mask shape (2008, 2000)\n",
      "FOV 021\n",
      "mask shape (2008, 2000)\n",
      "F021, no mask needed\n",
      "FOV 022\n",
      "mask shape (2008, 2000)\n",
      "F022, no mask needed\n",
      "FOV 023\n",
      "mask shape (2008, 2000)\n",
      "F023, no mask needed\n",
      "FOV 024\n",
      "mask shape (2008, 2000)\n",
      "F024, no mask needed\n",
      "FOV 026\n",
      "mask shape (2008, 2000)\n",
      "F026, no mask needed\n",
      "FOV 027\n",
      "mask shape (2008, 2000)\n",
      "F027, no mask needed\n",
      "FOV 029\n",
      "mask shape (2008, 2000)\n",
      "F029, no mask needed\n",
      "FOV 030\n",
      "mask shape (2008, 2000)\n",
      "F030, no mask needed\n",
      "FOV 031\n",
      "mask shape (2008, 2000)\n",
      "F031, no mask needed\n",
      "FOV 032\n",
      "mask shape (2008, 2000)\n",
      "FOV 033\n",
      "mask shape (2008, 2000)\n",
      "FOV 034\n",
      "mask shape (2008, 2000)\n",
      "F034, no mask needed\n",
      "FOV 035\n",
      "mask shape (2008, 2000)\n",
      "F035, no mask needed\n",
      "FOV 036\n",
      "mask shape (2008, 2000)\n",
      "FOV 038\n",
      "mask shape (2008, 2000)\n",
      "F038, no mask needed\n",
      "FOV 039\n",
      "mask shape (2008, 2000)\n",
      "FOV 040\n",
      "mask shape (2008, 2000)\n",
      "FOV 041\n",
      "mask shape (2008, 2000)\n",
      "F041, no mask needed\n",
      "FOV 042\n",
      "mask shape (2008, 2000)\n",
      "F042, no mask needed\n",
      "FOV 043\n",
      "mask shape (2008, 2000)\n",
      "F043, no mask needed\n",
      "FOV 044\n",
      "mask shape (2008, 2000)\n",
      "F044, no mask needed\n",
      "FOV 045\n",
      "mask shape (2008, 2000)\n",
      "F045, no mask needed\n",
      "FOV 046\n",
      "mask shape (2008, 2000)\n",
      "FOV 047\n",
      "mask shape (2008, 2000)\n",
      "F047, no mask needed\n",
      "FOV 048\n",
      "mask shape (2008, 2000)\n",
      "F048, no mask needed\n",
      "FOV 049\n",
      "mask shape (2008, 2000)\n",
      "FOV 051\n",
      "mask shape (2008, 2000)\n",
      "F051, no mask needed\n",
      "FOV 053\n",
      "mask shape (2008, 2000)\n",
      "F053, no mask needed\n",
      "FOV 054\n",
      "mask shape (2008, 2000)\n",
      "F054, no mask needed\n",
      "FOV 055\n",
      "mask shape (2008, 2000)\n",
      "F055, no mask needed\n",
      "FOV 056\n",
      "mask shape (2008, 2000)\n",
      "F056, no mask needed\n",
      "FOV 057\n",
      "mask shape (2008, 2000)\n",
      "F057, no mask needed\n",
      "FOV 058\n",
      "mask shape (2008, 2000)\n",
      "F058, no mask needed\n",
      "FOV 059\n",
      "mask shape (2008, 2000)\n",
      "F059, no mask needed\n",
      "FOV 060\n",
      "mask shape (2008, 2000)\n",
      "FOV 061\n",
      "mask shape (2008, 2000)\n",
      "F061, no mask needed\n",
      "FOV 063\n",
      "mask shape (2008, 2000)\n",
      "F063, no mask needed\n",
      "FOV 064\n",
      "mask shape (2008, 2000)\n",
      "F064, no mask needed\n",
      "FOV 065\n",
      "mask shape (2008, 2000)\n",
      "F065, no mask needed\n",
      "FOV 066\n",
      "mask shape (2008, 2000)\n",
      "F066, no mask needed\n",
      "FOV 067\n",
      "mask shape (2008, 2000)\n",
      "F067, no mask needed\n",
      "FOV 068\n",
      "mask shape (2008, 2000)\n",
      "F068, no mask needed\n",
      "FOV 069\n",
      "mask shape (2008, 2000)\n",
      "F069, no mask needed\n",
      "FOV 070\n",
      "mask shape (2008, 2000)\n",
      "F070, no mask needed\n",
      "FOV 071\n",
      "mask shape (2008, 2000)\n",
      "F071, no mask needed\n",
      "FOV 072\n",
      "mask shape (2008, 2000)\n",
      "F072, no mask needed\n",
      "FOV 073\n",
      "mask shape (2008, 2000)\n",
      "FOV 074\n",
      "mask shape (2008, 2000)\n",
      "F074, no mask needed\n",
      "FOV 075\n",
      "mask shape (2008, 2000)\n",
      "F075, no mask needed\n",
      "FOV 076\n",
      "mask shape (2008, 2000)\n",
      "F076, no mask needed\n",
      "FOV 078\n",
      "mask shape (2008, 2000)\n",
      "F078, no mask needed\n",
      "FOV 079\n",
      "mask shape (2008, 2000)\n",
      "F079, no mask needed\n",
      "FOV 080\n",
      "mask shape (2008, 2000)\n",
      "F080, no mask needed\n",
      "FOV 081\n",
      "mask shape (2008, 2000)\n",
      "F081, no mask needed\n",
      "FOV 083\n",
      "mask shape (2008, 2000)\n",
      "F083, no mask needed\n",
      "FOV 084\n",
      "mask shape (2008, 2000)\n",
      "F084, no mask needed\n",
      "FOV 085\n",
      "mask shape (2008, 2000)\n",
      "F085, no mask needed\n",
      "FOV 086\n",
      "mask shape (2008, 2000)\n",
      "FOV 087\n",
      "mask shape (2008, 2000)\n",
      "F087, no mask needed\n",
      "FOV 088\n",
      "mask shape (2008, 2000)\n",
      "F088, no mask needed\n",
      "FOV 089\n",
      "mask shape (2008, 2000)\n",
      "FOV 090\n",
      "mask shape (2008, 2000)\n",
      "F090, no mask needed\n",
      "FOV 091\n",
      "mask shape (2008, 2000)\n",
      "F091, no mask needed\n",
      "FOV 092\n",
      "mask shape (2008, 2000)\n",
      "F092, no mask needed\n",
      "FOV 093\n",
      "mask shape (2008, 2000)\n",
      "F093, no mask needed\n",
      "FOV 094\n",
      "mask shape (2008, 2000)\n",
      "F094, no mask needed\n",
      "FOV 095\n",
      "mask shape (2008, 2000)\n",
      "F095, no mask needed\n",
      "FOV 096\n",
      "mask shape (2008, 2000)\n",
      "F096, no mask needed\n",
      "FOV 097\n",
      "mask shape (2008, 2000)\n",
      "F097, no mask needed\n",
      "FOV 100\n",
      "mask shape (2008, 2000)\n",
      "F100, no mask needed\n",
      "FOV 101\n",
      "mask shape (2008, 2000)\n",
      "F101, no mask needed\n",
      "FOV 102\n",
      "mask shape (2008, 2000)\n",
      "F102, no mask needed\n",
      "FOV 103\n",
      "mask shape (2008, 2000)\n",
      "FOV 104\n",
      "mask shape (2008, 2000)\n",
      "F104, no mask needed\n",
      "FOV 105\n",
      "mask shape (2008, 2000)\n",
      "F105, no mask needed\n",
      "FOV 106\n",
      "mask shape (2008, 2000)\n",
      "F106, no mask needed\n",
      "FOV 107\n",
      "mask shape (2008, 2000)\n",
      "F107, no mask needed\n",
      "FOV 108\n",
      "mask shape (2008, 2000)\n",
      "FOV 109\n",
      "mask shape (2008, 2000)\n",
      "F109, no mask needed\n",
      "FOV 110\n",
      "mask shape (2008, 2000)\n",
      "F110, no mask needed\n",
      "FOV 111\n",
      "mask shape (2008, 2000)\n",
      "F111, no mask needed\n",
      "FOV 112\n",
      "mask shape (2008, 2000)\n",
      "F112, no mask needed\n",
      "FOV 113\n",
      "mask shape (2008, 2000)\n",
      "FOV 114\n",
      "mask shape (2008, 2000)\n",
      "F114, no mask needed\n",
      "FOV 115\n",
      "mask shape (2008, 2000)\n",
      "F115, no mask needed\n",
      "FOV 116\n",
      "mask shape (2008, 2000)\n",
      "F116, no mask needed\n",
      "FOV 117\n",
      "mask shape (2008, 2000)\n",
      "FOV 118\n",
      "mask shape (2008, 2000)\n",
      "F118, no mask needed\n",
      "FOV 119\n",
      "mask shape (2008, 2000)\n",
      "F119, no mask needed\n",
      "FOV 120\n",
      "mask shape (2008, 2000)\n",
      "F120, no mask needed\n",
      "FOV 121\n",
      "mask shape (2008, 2000)\n",
      "F121, no mask needed\n",
      "FOV 122\n",
      "mask shape (2008, 2000)\n",
      "F122, no mask needed\n",
      "FOV 123\n",
      "mask shape (2008, 2000)\n",
      "F123, no mask needed\n",
      "FOV 124\n",
      "mask shape (2008, 2000)\n",
      "F124, no mask needed\n",
      "FOV 125\n",
      "mask shape (2008, 2000)\n",
      "FOV 126\n",
      "mask shape (2008, 2000)\n",
      "F126, no mask needed\n",
      "FOV 127\n",
      "mask shape (2008, 2000)\n",
      "F127, no mask needed\n",
      "FOV 128\n",
      "mask shape (2008, 2000)\n",
      "F128, no mask needed\n",
      "FOV 129\n",
      "mask shape (2008, 2000)\n",
      "F129, no mask needed\n",
      "FOV 130\n",
      "mask shape (2008, 2000)\n",
      "F130, no mask needed\n",
      "FOV 131\n",
      "mask shape (2008, 2000)\n",
      "F131, no mask needed\n",
      "FOV 132\n",
      "mask shape (2008, 2000)\n",
      "F132, no mask needed\n",
      "FOV 133\n",
      "mask shape (2008, 2000)\n",
      "F133, no mask needed\n",
      "FOV 134\n",
      "mask shape (2008, 2000)\n",
      "F134, no mask needed\n",
      "FOV 135\n",
      "mask shape (2008, 2000)\n",
      "F135, no mask needed\n",
      "FOV 136\n",
      "mask shape (2008, 2000)\n",
      "F136, no mask needed\n",
      "FOV 137\n",
      "mask shape (2008, 2000)\n",
      "F137, no mask needed\n",
      "FOV 138\n",
      "mask shape (2008, 2000)\n",
      "FOV 139\n",
      "mask shape (2008, 2000)\n",
      "FOV 140\n",
      "mask shape (2008, 2000)\n",
      "F140, no mask needed\n",
      "FOV 141\n",
      "mask shape (2008, 2000)\n",
      "F141, no mask needed\n",
      "FOV 142\n",
      "mask shape (2008, 2000)\n",
      "F142, no mask needed\n",
      "FOV 143\n",
      "mask shape (2008, 2000)\n",
      "F143, no mask needed\n",
      "FOV 144\n",
      "mask shape (2008, 2000)\n",
      "F144, no mask needed\n",
      "FOV 145\n",
      "mask shape (2008, 2000)\n",
      "F145, no mask needed\n",
      "FOV 146\n",
      "mask shape (2008, 2000)\n",
      "F146, no mask needed\n",
      "FOV 147\n",
      "mask shape (2008, 2000)\n",
      "F147, no mask needed\n",
      "FOV 148\n",
      "mask shape (2008, 2000)\n",
      "F148, no mask needed\n",
      "FOV 149\n",
      "mask shape (2008, 2000)\n",
      "F149, no mask needed\n",
      "FOV 150\n",
      "mask shape (2008, 2000)\n",
      "FOV 151\n",
      "mask shape (2008, 2000)\n",
      "F151, no mask needed\n",
      "FOV 152\n",
      "mask shape (2008, 2000)\n",
      "F152, no mask needed\n",
      "FOV 153\n",
      "mask shape (2008, 2000)\n",
      "F153, no mask needed\n",
      "FOV 154\n",
      "mask shape (2008, 2000)\n",
      "FOV 155\n",
      "mask shape (2008, 2000)\n",
      "F155, no mask needed\n",
      "FOV 156\n",
      "mask shape (2008, 2000)\n",
      "F156, no mask needed\n",
      "FOV 158\n",
      "mask shape (2008, 2000)\n",
      "F158, no mask needed\n",
      "FOV 159\n",
      "mask shape (2008, 2000)\n",
      "F159, no mask needed\n",
      "FOV 160\n",
      "mask shape (2008, 2000)\n",
      "F160, no mask needed\n",
      "FOV 161\n",
      "mask shape (2008, 2000)\n",
      "F161, no mask needed\n",
      "FOV 162\n",
      "mask shape (2008, 2000)\n",
      "F162, no mask needed\n",
      "FOV 163\n",
      "mask shape (2008, 2000)\n",
      "F163, no mask needed\n",
      "FOV 164\n",
      "mask shape (2008, 2000)\n",
      "F164, no mask needed\n",
      "FOV 165\n",
      "mask shape (2008, 2000)\n",
      "FOV 166\n",
      "mask shape (2008, 2000)\n",
      "F166, no mask needed\n",
      "FOV 167\n",
      "mask shape (2008, 2000)\n",
      "F167, no mask needed\n",
      "FOV 168\n",
      "mask shape (2008, 2000)\n",
      "F168, no mask needed\n",
      "FOV 169\n",
      "mask shape (2008, 2000)\n",
      "F169, no mask needed\n",
      "FOV 170\n",
      "mask shape (2008, 2000)\n",
      "F170, no mask needed\n",
      "FOV 171\n",
      "mask shape (2008, 2000)\n",
      "F171, no mask needed\n",
      "FOV 172\n",
      "mask shape (2008, 2000)\n",
      "F172, no mask needed\n",
      "FOV 173\n",
      "mask shape (2008, 2000)\n",
      "F173, no mask needed\n",
      "FOV 174\n",
      "mask shape (2008, 2000)\n",
      "FOV 175\n",
      "mask shape (2008, 2000)\n",
      "F175, no mask needed\n",
      "FOV 176\n",
      "mask shape (2008, 2000)\n",
      "F176, no mask needed\n",
      "FOV 177\n",
      "mask shape (2008, 2000)\n",
      "F177, no mask needed\n",
      "FOV 178\n",
      "mask shape (2008, 2000)\n",
      "F178, no mask needed\n",
      "FOV 179\n",
      "mask shape (2008, 2000)\n",
      "F179, no mask needed\n",
      "FOV 180\n",
      "mask shape (2008, 2000)\n",
      "F180, no mask needed\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FOV 181\n",
      "mask shape (2008, 2000)\n",
      "FOV 182\n",
      "mask shape (2008, 2000)\n",
      "F182, no mask needed\n",
      "FOV 183\n",
      "mask shape (2008, 2000)\n",
      "F183, no mask needed\n",
      "FOV 184\n",
      "mask shape (2008, 2000)\n",
      "FOV 185\n",
      "mask shape (2008, 2000)\n",
      "F185, no mask needed\n",
      "FOV 186\n",
      "mask shape (2008, 2000)\n",
      "FOV 188\n",
      "mask shape (2008, 2000)\n",
      "F188, no mask needed\n",
      "FOV 189\n",
      "mask shape (2008, 2000)\n",
      "F189, no mask needed\n",
      "FOV 190\n",
      "mask shape (2008, 2000)\n",
      "F190, no mask needed\n",
      "FOV 191\n",
      "mask shape (2008, 2000)\n",
      "F191, no mask needed\n",
      "FOV 192\n",
      "mask shape (2008, 2000)\n",
      "F192, no mask needed\n",
      "FOV 193\n",
      "mask shape (2008, 2000)\n",
      "F193, no mask needed\n",
      "FOV 194\n",
      "mask shape (2008, 2000)\n",
      "F194, no mask needed\n",
      "FOV 195\n",
      "mask shape (2008, 2000)\n",
      "F195, no mask needed\n",
      "FOV 196\n",
      "mask shape (2008, 2000)\n",
      "F196, no mask needed\n",
      "FOV 197\n",
      "mask shape (2008, 2000)\n",
      "FOV 199\n",
      "mask shape (2008, 2000)\n",
      "F199, no mask needed\n",
      "FOV 201\n",
      "mask shape (2008, 2000)\n",
      "F201, no mask needed\n",
      "FOV 202\n",
      "mask shape (2008, 2000)\n",
      "F202, no mask needed\n",
      "FOV 207\n",
      "mask shape (2008, 2000)\n",
      "F207, no mask needed\n",
      "FOV 208\n",
      "mask shape (2008, 2000)\n",
      "F208, no mask needed\n",
      "FOV 209\n",
      "mask shape (2008, 2000)\n",
      "FOV 210\n",
      "mask shape (2008, 2000)\n",
      "F210, no mask needed\n",
      "FOV 211\n",
      "mask shape (2008, 2000)\n",
      "F211, no mask needed\n",
      "FOV 212\n",
      "mask shape (2008, 2000)\n",
      "F212, no mask needed\n",
      "FOV 217\n",
      "mask shape (2008, 2000)\n",
      "F217, no mask needed\n",
      "FOV 218\n",
      "mask shape (2008, 2000)\n",
      "F218, no mask needed\n",
      "FOV 219\n",
      "mask shape (2008, 2000)\n",
      "F219, no mask needed\n",
      "FOV 220\n",
      "mask shape (2008, 2000)\n",
      "F220, no mask needed\n",
      "FOV 221\n",
      "mask shape (2008, 2000)\n",
      "F221, no mask needed\n"
     ]
    }
   ],
   "source": [
    "NUM_FOVS = 195\n",
    "merged = iter(glob.glob('merged/*')) \n",
    "for fov in range(NUM_FOVS):\n",
    "    merged_name = next(merged)\n",
    "    img = imread(merged_name)\n",
    "    img = img.astype(np.uint16)\n",
    "    FOV_num = merged_name.split('/F')[1][0:3]\n",
    "    print(\"FOV\", FOV_num)\n",
    "    \n",
    "    img = img.transpose(1,0,2,3) \n",
    "    # C Z Y X format. you may use image as it is but should change bottom line to axis=(0,2,3)\n",
    "    count0 = np.count_nonzero(img==0, axis=(1,2,3))\n",
    "    count20 = np.count_nonzero(img<20, axis=(1,2,3))\n",
    "    count65K = np.count_nonzero(img>65000, axis=(1,2,3))\n",
    "\n",
    "    # initialize mask, save_im_flag\n",
    "    save_im_flag = False\n",
    "    mask = np.zeros(img[0,0,...].shape)\n",
    "    print(\"mask shape\", mask.shape)\n",
    "    mask = mask.astype('bool')\n",
    "    # check for each channel;\n",
    "    for ch in range(count0.shape[0]):\n",
    "        # record in dataframe\n",
    "        count0_df.loc[FOV_num, ch] = count0[ch]\n",
    "        count20_df.loc[FOV_num, ch] = count20[ch]\n",
    "        count65K_df.loc[FOV_num, ch] = count65K[ch]\n",
    "        if ch in final_inds:\n",
    "            # whether to make mask;\n",
    "            if count20[ch] > 2000: # if there are more than 2000 <20 values\n",
    "                save_im_flag = True\n",
    "                for z in range(img.shape[1]): # second index is Z dimension, after Transposing\n",
    "                    mask = mask | (img[ch, z] < 20).astype('bool')\n",
    "            exp_mask = expand_labels(mask, distance=5) # expand mask with 5 pixel padding\n",
    "\n",
    "            # check if area is large enough? 200\n",
    "            labels = measure.label(mask)\n",
    "            _df = measure.regionprops_table(labels, mask, properties=['label','area','centroid','bbox'])\n",
    "            df=pd.DataFrame(_df)\n",
    "            df=df.set_index('label')\n",
    "            if df.area.max() < 200: # if area is less than 200, turn off flag and don't do anything?\n",
    "                save_im_flag = False\n",
    "\n",
    "\n",
    "    # export mask if flag is true\n",
    "    if save_im_flag:\n",
    "        #sFOV = str(fov).zfill(3)\n",
    "        fname = f'{OUT_DIR}/F{FOV_num}_mask.tif'\n",
    "        tifffile.imwrite(fname, exp_mask.astype('uint8'), imagej = True, photometric='minisblack', metadata={'axes':'YX'})\n",
    "    else:\n",
    "        print(f'F{FOV_num}, no mask needed') # export later!\n",
    "        \n",
    "    \n",
    "# save dataframe \n",
    "count0_df.to_csv('count0_df.csv', sep=',')\n",
    "count20_df.to_csv('count20_df.csv', sep=',')\n",
    "count65K_df.to_csv('count65K_df.csv', sep=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mask percentage >3% and <5%\n",
      "014\n",
      "3.968874501992032\n",
      "mask percentage >3% and <5%\n",
      "033\n",
      "3.3864541832669324\n",
      "mask percentage >5% and <7%\n",
      "049\n",
      "6.954631474103586\n",
      "mask percentage >3% and <5%\n",
      "060\n",
      "4.058217131474104\n",
      "mask percentage >3% and <5%\n",
      "108\n",
      "3.592106573705179\n",
      "mask percentage >3% and <5%\n",
      "117\n",
      "3.744546812749004\n",
      "mask percentage >3% and <5%\n",
      "125\n",
      "4.677739043824701\n",
      "mask percentage >3% and <5%\n",
      "138\n",
      "3.127241035856574\n",
      "mask percentage >3% and <5%\n",
      "181\n",
      "4.229506972111554\n",
      "mask percentage >3% and <5%\n",
      "186\n",
      "4.780901394422311\n",
      "mask percentage >3% and <5%\n",
      "197\n",
      "4.3773406374501995\n",
      "mask percentage >3% and <5%\n",
      "209\n",
      "4.443227091633466\n"
     ]
    }
   ],
   "source": [
    "NUM_MASKS = 31\n",
    "masks = iter(glob.glob('mask/*')) \n",
    "for fov in range(NUM_MASKS):\n",
    "    mask_name = next(masks)\n",
    "    img = imread(mask_name)\n",
    "    img = img.astype(np.uint16)\n",
    "    FOV_num = mask_name.split('/F')[1][0:3]\n",
    "    \n",
    "    mask_perc = 100*np.count_nonzero(img == 1)/(img.shape[0]*img.shape[1])\n",
    "    \n",
    "    if mask_perc > 7:\n",
    "        print(\"mask percentage >7%\")\n",
    "        print(FOV_num)\n",
    "        print(mask_perc)\n",
    "    if mask_perc > 5 and mask_perc < 7:\n",
    "        print(\"mask percentage >5% and <7%\")\n",
    "        print(FOV_num)\n",
    "        print(mask_perc)\n",
    "    if mask_perc > 3 and mask_perc < 5:\n",
    "        print(\"mask percentage >3% and <5%\")\n",
    "        print(FOV_num)\n",
    "        print(mask_perc)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Max Projections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as patches\n",
    "from pathlib import Path\n",
    "import numpy as np\n",
    "import os\n",
    "import sys\n",
    "import glob\n",
    "from imageio import volread as imread"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mkdir max"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FOV 001\n",
      "Prepare CP inputs - max projection only\n",
      "Saving MAX... ./max/fname\n",
      "(25, 2008, 2000)\n",
      "FOV 002\n",
      "Prepare CP inputs - max projection only\n",
      "Saving MAX... ./max/fname\n",
      "(25, 2008, 2000)\n",
      "FOV 004\n",
      "Prepare CP inputs - max projection only\n",
      "Saving MAX... ./max/fname\n",
      "(25, 2008, 2000)\n",
      "FOV 005\n",
      "Prepare CP inputs - max projection only\n",
      "Saving MAX... ./max/fname\n",
      "(25, 2008, 2000)\n",
      "FOV 006\n",
      "Prepare CP inputs - max projection only\n",
      "Saving MAX... ./max/fname\n",
      "(25, 2008, 2000)\n",
      "FOV 007\n",
      "Prepare CP inputs - max projection only\n",
      "Saving MAX... ./max/fname\n",
      "(25, 2008, 2000)\n",
      "FOV 008\n",
      "Prepare CP inputs - max projection only\n",
      "Saving MAX... ./max/fname\n",
      "(25, 2008, 2000)\n",
      "FOV 009\n",
      "Prepare CP inputs - max projection only\n",
      "Saving MAX... ./max/fname\n",
      "(25, 2008, 2000)\n",
      "FOV 010\n",
      "Prepare CP inputs - max projection only\n",
      "Saving MAX... ./max/fname\n",
      "(25, 2008, 2000)\n",
      "FOV 012\n",
      "Prepare CP inputs - max projection only\n",
      "Saving MAX... ./max/fname\n",
      "(25, 2008, 2000)\n",
      "FOV 014\n",
      "Prepare CP inputs - max projection only\n",
      "Saving MAX... ./max/fname\n",
      "(25, 2008, 2000)\n",
      "FOV 015\n",
      "Prepare CP inputs - max projection only\n",
      "Saving MAX... ./max/fname\n",
      "(25, 2008, 2000)\n",
      "FOV 016\n",
      "Prepare CP inputs - max projection only\n",
      "Saving MAX... ./max/fname\n",
      "(25, 2008, 2000)\n",
      "FOV 017\n",
      "Prepare CP inputs - max projection only\n",
      "Saving MAX... ./max/fname\n",
      "(25, 2008, 2000)\n",
      "FOV 019\n",
      "Prepare CP inputs - max projection only\n",
      "Saving MAX... ./max/fname\n",
      "(25, 2008, 2000)\n"
     ]
    }
   ],
   "source": [
    "MAX_DIR = 'max'\n",
    "merged = iter(glob.glob('merged/*')) \n",
    "for FOV in range(NUM_FOVS):\n",
    "    merged_name = next(merged)\n",
    "    img = imread(merged_name)\n",
    "    img = img.astype(np.uint16)\n",
    "    FOV_num = merged_name.split('/F')[1][0:3]\n",
    "    print(\"FOV\", FOV_num)\n",
    "    \n",
    "    print(\"Prepare CP inputs - max projection only\")\n",
    "    img_max = img.max(0) # take max projection # ZCYX\n",
    "    im_to_save = img_max[final_inds].copy()\n",
    "    fname = f'F{FOV_num}_max.tif'\n",
    "    \n",
    "    print('Saving MAX...', f'./{MAX_DIR}/fname') \n",
    "    print(im_to_save.shape)\n",
    "    tifffile.imwrite(f'./{MAX_DIR}/'+fname, im_to_save, imagej = True,\n",
    "                    photometric='minisblack', metadata={'axes':'CYX'})"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
