20220103_IMC_pipeline.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78a56982",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import libraries\n",
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import shutil\n",
    "import matplotlib.pyplot as plt\n",
    "import re\n",
    "from skimage import io\n",
    "import tifffile\n",
    "from scipy.ndimage import median_filter\n",
    "\n",
    "# Set Paths\n",
    "#codedir = '/home/groups/graylab_share/OMERO.rdsStore/engje/Data/cmIF'\n",
    "codedir = os.getcwd()\n",
    "rootdir = f'{codedir}'\n",
    "#tiffdir = f'{rootdir}/RawImages'\n",
    "qcdir = f'{rootdir}/QC'\n",
    "#regdir = f'{rootdir}/RegisteredImages'\n",
    "subdir = f'{rootdir}/SubtractedRegisteredImages'\n",
    "segdir = f'{rootdir}/Segmentation'\n",
    "#czidir = f'{rootdir}/czis'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18ba42c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start Preprocessing\n",
    "os.chdir('/home/groups/graylab_share/data/engje/Data/')\n",
    "from mplex_image import preprocess, mpimage, cmif\n",
    "os.chdir(codedir)\n",
    "preprocess.cmif_mkdir([qcdir,segdir,subdir])\n",
    "\n",
    "# List slides to be processed\n",
    "ls_sample = [\n",
    " ]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f1ef3f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy\n",
    "from skimage.filters import unsharp_mask\n",
    "import skimage\n",
    "from skimage.restoration import (denoise_tv_chambolle, denoise_bilateral,\n",
    "                                 denoise_wavelet, estimate_sigma)\n",
    "from skimage import color, morphology\n",
    "from skimage.transform import rescale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "693da2db",
   "metadata": {},
   "outputs": [],
   "source": [
    "codedir"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "221db693",
   "metadata": {},
   "source": [
    "# image renaming"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d81d6a8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(f'{codedir}/ome/ome/')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17b2ffc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_s = pd.read_csv(f'{codedir}/Basel_Zuri_StainingPanel.csv')\n",
    "ls_drop = [27,42,47,49]\n",
    "df_s.drop(ls_drop,inplace=True)\n",
    "#naming errors\n",
    "df_s.loc[37,'Target'] = 'Ecad'\n",
    "df_s.loc[26,'Target'] = 'ER'\n",
    "ls_marker = ['0ArArArAr80Di', 'utheni_10311239Ru96Di', 'utheni_10311240Ru98Di', 'utheni_10311241Ru99Di', 'utheni_10311242Ru100Di',\n",
    " 'utheni_10311243Ru101Di', 'utheni_10311244Ru102Di', 'utheni_10311245Ru104Di', 'istone_1261726In113Di', '15InIn115Di',\n",
    " '134XeXe134Di', 'Histone_473968La139Di', 'Cytoker_651779Pr141Di', 'Fibrone_3281668Nd142Di', 'Cytoker_3111576Nd143Di',\n",
    " 'Cytoker_971099Nd144Di', 'UnknownProteinNd145Di', 'CD68_77877Nd146Di', 'Keratin_346876Sm147Di', 'SMA_174864Nd148Di',\n",
    " 'Vimenti_1921755Sm149Di', 'c-Myc_322787Nd150Di', 'c-erbB-_201487Eu151Di', 'CD3epsi_8001752Sm152Di', 'Histone_phospho_383738Eu153Di',\n",
    " 'Erk12_phospho_105894Sm154Di', 'Slug_3521227Gd155Di', 'Estroge_112475Gd156Di', 'Progest_312878Gd158Di', 'p53_207736Tb159Di',\n",
    " 'CD44_6967Gd160Di', 'EpCAM_2971330Dy161Di', 'CD45_71790Dy162Di', 'GATA3_117792Dy163Di', 'CD20_361077Dy164Di',\n",
    " 'b-Caten_1971527Ho165Di', 'Carboni_92964Er166Di', 'E-Cadhe_1031747Er167Di', 'Ki-67_1441101Er168Di', 'EGFR_1021522Tm169Di',\n",
    " 'S6_phospho_166866Er170Di', 'Sox9_483739Yb171Di', 'vWF_378871Yb172Di', 'mTOR_phospho_204417Yb173Di', 'Cytoker_98922Yb174Di',\n",
    " 'panCyto_234832Lu175Di', 'cleaved_198883Yb176Di', 'Iridium_10331253Ir191Di', 'Iridium_10331254Ir193Di', '202HgHg202Di',\n",
    " '204PbPb204Di', '206PbPb206Di']\n",
    "#ls_metal = [item.replace('_phospho','-phospho').replace('Ru9','Ru09').replace('Ar80','Ar080').split('Di')[0][-5::] for item in ls_marker]\n",
    "ls_metal = [item.split('Di')[0][-5::].replace('rAr80','Ar80').replace('9Ru','Ru').replace('0Ru','Ru').replace('1Ru','Ru')  for item in ls_marker]\n",
    "d_metal = dict(zip(ls_marker,ls_metal))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fad57071",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#ER is Rbt IgG :/; no its not!\n",
    "#Registered-R0_R0c2.R0c3.R0c4.R0c5_BC44290-146-Scene-01_c1_ORG.tif\n",
    "df_result = pd.DataFrame()\n",
    "se_cyto = pd.Series(['Ecad','Vimentin','CD44','CD45'])\n",
    "#uint16 0 to 65535\n",
    "for s_file in os.listdir():\n",
    "    if s_file.find('full.tiff') > -1:\n",
    "        if s_file.find('BaselTMA') > -1:\n",
    "            print(s_file)\n",
    "            s_slide = s_file.split('_')[7]\n",
    "            s_scene = s_file.split('_')[8]\n",
    "            with tifffile.TiffFile(s_file) as tif:\n",
    "                array = tif.asarray()\n",
    "                omexml_string = tif.ome_metadata\n",
    "                ls_unk = []\n",
    "                imgs_cyto = []\n",
    "                for idx in range(array.shape[0]):\n",
    "                    #idx=37\n",
    "                    img = array[idx]\n",
    "                    #dtype\n",
    "                    img = np.clip(img, a_min=0, a_max=65535).astype('uint16')\n",
    "                    i_begin = omexml_string.find(f'ID=\"Channel:0:{idx}\" Name=\"')\n",
    "                    i_end = omexml_string[i_begin:].find('\" SamplesPerPixel')\n",
    "                    s_marker = omexml_string[i_begin + 24:i_begin + i_end]\n",
    "                    #print(idx)\n",
    "                    #print(s_marker.split('_')[0])\n",
    "                    s_metal = d_metal[s_marker]\n",
    "                    if len(df_s[df_s.loc[:,'Metal Tag'] == s_metal].Target) == 1:\n",
    "                        s_target = df_s[df_s.loc[:,'Metal Tag'] == s_metal].Target.iloc[0]\n",
    "                    else:\n",
    "                        #print(f'Error in target {s_metal}')\n",
    "                        #print(df_s[df_s.loc[:,'Metal Tag'] == s_metal].Target)\n",
    "                        s_target = s_metal\n",
    "                        ls_unk.append(s_marker)\n",
    "                    s_name = s_target.replace(' ','').replace('/','').replace('(KRT14)','').replace('-','')\n",
    "                    #print(s_name)\n",
    "                    if s_name =='DNA1':\n",
    "                        dapi_scale = skimage.exposure.rescale_intensity(img,in_range=(np.quantile(img,0.03),1.5*np.quantile(img,0.99999)))\n",
    "                        s_filename = f\"Registered-R1_{s_name}_{s_slide}-Scene-{s_scene}_c1_ORG.tif\"\n",
    "                        #print(s_filename)\n",
    "                    elif s_name =='DNA2':\n",
    "                        s_filename = f\"Registered-R2_{s_name}_{s_slide}-Scene-{s_scene}_c1_ORG.tif\"\n",
    "                        #print(s_filename)\n",
    "                    elif se_cyto.str.contains(s_name).sum() > 0:\n",
    "                        cyto_scale = skimage.exposure.rescale_intensity(img,in_range=(np.quantile(img,0.03),1.5*np.quantile(img,0.99999)))\n",
    "                        imgs_cyto.append(cyto_scale)\n",
    "                    if s_name.find('DNA') == -1:\n",
    "                        s_filename = f\"Registered-R{idx}_{s_name}_{s_slide}-Scene-{s_scene}_c2_ORG.tif\"\n",
    "                    #hot pixels fix\n",
    "                    blurred = median_filter(img, size=2)\n",
    "                    difference = img.astype('int16') - blurred.astype('int16')\n",
    "                    threshold = 10*np.std(difference)\n",
    "                    i_filter = sum(img > int(threshold)).sum()/img.size\n",
    "                    df_result.loc[s_filename,s_name] = i_filter\n",
    "                    img[img > int(threshold)] = blurred[img > int(threshold)]\n",
    "                    #save\n",
    "                    io.imsave(f\"{subdir}/{s_filename}\", img,check_contrast=False)\n",
    "                    #break\n",
    "                #break\n",
    "                #dapis\n",
    "                '''\n",
    "                dapi = np.clip((dapi_scale),a_min=0,a_max=65535)\n",
    "                dapi_gamma = skimage.exposure.adjust_gamma(dapi, gamma=0.4)\n",
    "                \n",
    "                ##pretty good for cellpose, not good enough though\n",
    "                ##noisy = added_gamma[200:450,300:550]\n",
    "                ##tv = denoise_tv_chambolle(noisy, weight=0.25)\n",
    "                ##image_rescaled = rescale(tv, 2, anti_aliasing=False)\n",
    "                ##io.imsave(f'{segdir}/test/Registered-R1_DNA1r_{s_slide}-Scene-{s_scene}_c1.png',(image_rescaled*255).astype('uint8'))\n",
    "                \n",
    "                #The mesmer job expects images of shape[(-1, 256, 256, 2)]. use cyto mip and dapi\n",
    "                mip = np.stack(imgs_cyto).max(axis=0)\n",
    "                cyto = np.clip((mip),a_min=0,a_max=65535)\n",
    "                cyto_gamma = skimage.exposure.adjust_gamma(cyto, gamma=0.4)\n",
    "                #de-noise        \n",
    "                tv = denoise_tv_chambolle(dapi_gamma, weight=0.1)\n",
    "                cyto_tv = denoise_tv_chambolle(cyto_gamma, weight=0.05)\n",
    "                zdh = np.dstack((tv,cyto_tv,np.zeros(tv.shape)))\n",
    "                #io.imsave(f'{segdir}/Registered-R1_DNA1.Cyto_{s_slide}-Scene-{s_scene}_c1c2.png',(zdh*255).astype('uint8'))\n",
    "                '''\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28960cba",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_result.to_csv('percent_pixels_filtered.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e0c01bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "i_max = np.quantile(img,.9998)\n",
    "io.imshow(img[0:100,0:100],vmin=0, vmax=i_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0041234",
   "metadata": {},
   "outputs": [],
   "source": [
    "io.imshow(blurred[0:100,0:100],vmin=0, vmax=i_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4890661e",
   "metadata": {},
   "outputs": [],
   "source": [
    "io.imshow(difference[0:100,0:100],cmap='bwr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b81e7c5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15d7ed8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#hot pixels begone\n",
    "s_file = 'BaselTMA_SP43_15.9kx19.3ky_10000x5500_10_20170928_86_111_X8Y6_105_a0_full.tiff'\n",
    "idx=28\n",
    "with tifffile.TiffFile(s_file) as tif:\n",
    "    array = tif.asarray()\n",
    "    img = array[idx]\n",
    "i_max = np.quantile(img,.9998)\n",
    "io.imshow(img,vmin=0, vmax=i_max)\n",
    "#hot pixels begone\n",
    "blurred = median_filter(img, size=2)\n",
    "difference = img.astype('int16') - blurred.astype('int16')\n",
    "threshold = 10*np.std(difference)\n",
    "img[img > int(threshold)] = 0\n",
    "fig, ax = plt.subplots()\n",
    "io.imshow(img,vmin=0, vmax=i_max,ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef12e93e",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#zdh = np.dstack((np.zeros(tv.shape),cyto_tv,tv))\n",
    "#io.imshow(zdh)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52bc04dd",
   "metadata": {},
   "source": [
    "## segmentation example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "661ccbe0",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Download multiplex data\n",
    "from deepcell.datasets import multiplex_tissue\n",
    "((X_train, y_train),(X_test, y_test)) = multiplex_tissue.load_data()\n",
    "from deepcell.utils.plot_utils import create_rgb_image\n",
    "\n",
    "# create rgb overlay of image data for visualization\n",
    "rgb_images = create_rgb_image(X_train, channel_colors=['green', 'blue'])\n",
    "\n",
    "# select index for displaying\n",
    "idx = 6 #3\n",
    "from matplotlib import pyplot as plt\n",
    "fig, ax = plt.subplots(1, 3, figsize=(15, 15))\n",
    "ax[0].imshow(X_train[idx, ..., 0])\n",
    "ax[1].imshow(X_train[idx, ..., 1])\n",
    "ax[2].imshow(rgb_images[idx, ...])\n",
    "\n",
    "ax[0].set_title('Nuclear channel')\n",
    "ax[1].set_title('Membrane channel')\n",
    "ax[2].set_title('Overlay')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40206847",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/cuda/extras/CUPTI #enter in terminal\n",
    "from deepcell.applications import Mesmer\n",
    "app = Mesmer()\n",
    "print('Training Resolution:', app.model_mpp, 'microns per pixel')\n",
    "segmentation_predictions = app.predict(X_train, image_mpp=0.5)\n",
    "# create overlay of predictions\n",
    "from deepcell.utils.plot_utils import make_outline_overlay\n",
    "overlay_data = make_outline_overlay(rgb_data=rgb_images, predictions=segmentation_predictions)\n",
    "fig, ax = plt.subplots(1, 2, figsize=(15, 15))\n",
    "ax[0].imshow(rgb_images[idx, ...])\n",
    "ax[1].imshow(overlay_data[idx, ...])\n",
    "ax[0].set_title('Raw data')\n",
    "ax[1].set_title('Predictions')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4ccf77d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "segmentation_predictions_nuc = app.predict(X_train, image_mpp=0.5, compartment='nuclear')\n",
    "overlay_data_nuc = make_outline_overlay(rgb_data=rgb_images, predictions=segmentation_predictions_nuc)\n",
    "fig, ax = plt.subplots(1, 2, figsize=(15, 15))\n",
    "ax[0].imshow(rgb_images[idx, ...])\n",
    "ax[1].imshow(overlay_data_nuc[idx, ...])\n",
    "ax[0].set_title('Raw data')\n",
    "ax[1].set_title('Nuclear Predictions')\n",
    "plt.show()\n",
    "#rgb_images.shape\n",
    "#X_train.shape\n",
    "#rgb_images.dtype\n",
    "#rgb_images.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac0cf44f",
   "metadata": {},
   "source": [
    "# segment IMC data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "594e147d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import io, measure, segmentation\n",
    "from deepcell.applications import Mesmer\n",
    "from deepcell.utils.plot_utils import create_rgb_image\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dbfe2cfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_file = [#'Registered-R1_DNA1.Cyto_46-Scene-X10Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_43-Scene-X5Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_22-Scene-X12Y6_c1c2.png','Registered-R1_DNA1.Cyto_95-Scene-X14Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_213-Scene-X7Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_233-Scene-X13Y7_c1c2.png','Registered-R1_DNA1.Cyto_143-Scene-X1Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_252-Scene-X6Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_3-Scene-X13Y6_c1c2.png','Registered-R1_DNA1.Cyto_39-Scene-X1Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_239-Scene-X11Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_31-Scene-X15Y2_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_113-Scene-X2Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_80-Scene-X4Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_14-Scene-X1Y8_c1c2.png','Registered-R1_DNA1.Cyto_262-Scene-X10Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_229-Scene-X1Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_17-Scene-X12Y4_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_183-Scene-X3Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_83-Scene-X13Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_282-Scene-X2Y3_c1c2.png','Registered-R1_DNA1.Cyto_68-Scene-X9Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_96-Scene-X9Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_117-Scene-X13Y3_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_57-Scene-X8Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_238-Scene-X7Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_108-Scene-X13Y8_c1c2.png','Registered-R1_DNA1.Cyto_68-Scene-X10Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_111-Scene-X8Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_80-Scene-X3Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_231-Scene-X6Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_222-Scene-X7Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_174-Scene-X9Y6_c1c2.png','Registered-R1_DNA1.Cyto_8-Scene-X5Y9_c1c2.png', 'Registered-R1_DNA1.Cyto_242-Scene-X5Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_20-Scene-X12Y5_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_154-Scene-X8Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_137-Scene-X9Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_180-Scene-X6Y2_c1c2.png','Registered-R1_DNA1.Cyto_235-Scene-X11Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_133-Scene-X3Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_99-Scene-X5Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_61-Scene-X4Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_47-Scene-X16Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_143-Scene-X2Y2_c1c2.png','Registered-R1_DNA1.Cyto_156-Scene-X16Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_169-Scene-X3Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_49-Scene-X1Y5_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_141-Scene-X11Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_11-Scene-X13Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_215-Scene-X4Y1_c1c2.png','Registered-R1_DNA1.Cyto_131-Scene-X13Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_164-Scene-X9Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_117-Scene-X14Y3_c1c2.png','Registered-R1_DNA1.Cyto_136-Scene-X9Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_34-Scene-X12Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_250-Scene-X2Y3_c1c2.png','Registered-R1_DNA1.Cyto_69-Scene-X5Y3_c1c2.png','Registered-R1_DNA1.Cyto_239-Scene-X12Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_243-Scene-X14Y3_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_220-Scene-X10Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_162-Scene-X16Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_32-Scene-X4Y6_c1c2.png','Registered-R1_DNA1.Cyto_238-Scene-X8Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_91-Scene-X12Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_98-Scene-X11Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_185-Scene-X7Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_172-Scene-X9Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_224-Scene-X3Y7_c1c2.png','Registered-R1_DNA1.Cyto_95-Scene-X13Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_186-Scene-X6Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_158-Scene-X12Y6_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_173-Scene-X2Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_44-Scene-X2Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_45-Scene-X3Y3_c1c2.png','Registered-R1_DNA1.Cyto_267-Scene-X13Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_225-Scene-X2Y9_c1c2.png', 'Registered-R1_DNA1.Cyto_92-Scene-X12Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_40-Scene-X16Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_187-Scene-X3Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_25-Scene-X4Y2_c1c2.png','Registered-R1_DNA1.Cyto_120-Scene-X6Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_67-Scene-X14Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_250-Scene-X1Y3_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_196-Scene-X5Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_211-Scene-X4Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_74-Scene-X15Y3_c1c2.png','Registered-R1_DNA1.Cyto_105-Scene-X9Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_264-Scene-X14Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_64-Scene-X14Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_29-Scene-X7Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_64-Scene-X13Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_61-Scene-X3Y4_c1c2.png','Registered-R1_DNA1.Cyto_116-Scene-X4Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_152-Scene-X15Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_151-Scene-X14Y7_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_119-Scene-X9Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_104-Scene-X7Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_103-Scene-X6Y4_c1c2.png','Registered-R1_DNA1.Cyto_136-Scene-X10Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_274-Scene-X13Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_275-Scene-X5Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_209-Scene-X11Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_23-Scene-X15Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_84-Scene-X8Y1_c1c2.png','Registered-R1_DNA1.Cyto_54-Scene-X12Y5_c1c2.png','Registered-R1_DNA1.Cyto_77-Scene-X3Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_167-Scene-X14Y8_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_260-Scene-X9Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_81-Scene-X11Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_78-Scene-X3Y9_c1c2.png','Registered-R1_DNA1.Cyto_142-Scene-X5Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_205-Scene-X3Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_85-Scene-X9Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_274-Scene-X14Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_283-Scene-X2Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_90-Scene-X7Y6_c1c2.png','Registered-R1_DNA1.Cyto_218-Scene-X10Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_155-Scene-X4Y9_c1c2.png', 'Registered-R1_DNA1.Cyto_4-Scene-X15Y3_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_10-Scene-X1Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_132-Scene-X14Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_4-Scene-X16Y3_c1c2.png','Registered-R1_DNA1.Cyto_134-Scene-X3Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_160-Scene-X13Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_241-Scene-X8Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_56-Scene-X2Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_280-Scene-X9Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_246-Scene-X2Y6_c1c2.png','Registered-R1_DNA1.Cyto_70-Scene-X6Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_6-Scene-X12Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_271-Scene-X11Y8_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_263-Scene-X1Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_255-Scene-X12Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_278-Scene-X4Y2_c1c2.png','Registered-R1_DNA1.Cyto_209-Scene-X12Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_284-Scene-X8Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_16-Scene-X4Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_112-Scene-X5Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_186-Scene-X5Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_144-Scene-X16Y1_c1c2.png','Registered-R1_DNA1.Cyto_261-Scene-X11Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_172-Scene-X10Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_165-Scene-X5Y6_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_148-Scene-X5Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_19-Scene-X16Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_244-Scene-X3Y1_c1c2.png','Registered-R1_DNA1.Cyto_106-Scene-X2Y9_c1c2.png','Registered-R1_DNA1.Cyto_2-Scene-X2Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_21-Scene-X3Y5_c1c2.png','Registered-R1_DNA1.Cyto_141-Scene-X12Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_258-Scene-X4Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_17-Scene-X11Y4_c1c2.png','Registered-R1_DNA1.Cyto_215-Scene-X3Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_156-Scene-X15Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_228-Scene-X7Y3_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_214-Scene-X15Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_171-Scene-X2Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_5-Scene-X8Y7_c1c2.png','Registered-R1_DNA1.Cyto_37-Scene-X7Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_168-Scene-X10Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_59-Scene-X3Y9_c1c2.png', 'Registered-R1_DNA1.Cyto_270-Scene-X6Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_55-Scene-X10Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_96-Scene-X10Y3_c1c2.png','Registered-R1_DNA1.Cyto_33-Scene-X14Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_170-Scene-X3Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_19-Scene-X15Y2_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_254-Scene-X11Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_53-Scene-X6Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_145-Scene-X2Y4_c1c2.png','Registered-R1_DNA1.Cyto_197-Scene-X4Y9_c1c2.png', 'Registered-R1_DNA1.Cyto_142-Scene-X6Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_210-Scene-X4Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_122-Scene-X12Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_170-Scene-X4Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_72-Scene-X12Y1_c1c2.png','Registered-R1_DNA1.Cyto_192-Scene-X8Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_135-Scene-X8Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_127-Scene-X2Y8_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_244-Scene-X4Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_84-Scene-X7Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_130-Scene-X1Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_206-Scene-X1Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_245-Scene-X6Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_176-Scene-X3Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_79-Scene-X7Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_236-Scene-X11Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_62-Scene-X7Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_36-Scene-X5Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_277-Scene-X12Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_279-Scene-X14Y6_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_201-Scene-X14Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_104-Scene-X8Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_284-Scene-X7Y2_c1c2.png','Registered-R1_DNA1.Cyto_83-Scene-X14Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_237-Scene-X16Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_195-Scene-X11Y6_c1c2.png','Registered-R1_DNA1.Cyto_107-Scene-X4Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_230-Scene-X11Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_16-Scene-X3Y4_c1c2.png','Registered-R1_DNA1.Cyto_15-Scene-X1Y9_c1c2.png', 'Registered-R1_DNA1.Cyto_13-Scene-X5Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_285-Scene-X7Y7_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_259-Scene-X10Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_257-Scene-X3Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_199-Scene-X8Y5_c1c2.png','Registered-R1_DNA1.Cyto_137-Scene-X10Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_253-Scene-X14Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_50-Scene-X8Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_260-Scene-X10Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_159-Scene-X11Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_38-Scene-X4Y7_c1c2.png','Registered-R1_DNA1.Cyto_128-Scene-X16Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_82-Scene-X2Y9_c1c2.png', 'Registered-R1_DNA1.Cyto_88-Scene-X11Y5_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_278-Scene-X3Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_116-Scene-X3Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_91-Scene-X11Y3_c1c2.png','Registered-R1_DNA1.Cyto_269-Scene-X5Y9_c1c2.png', 'Registered-R1_DNA1.Cyto_144-Scene-X15Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_207-Scene-X1Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_105-Scene-X10Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_150-Scene-X6Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_42-Scene-X14Y5_c1c2.png','Registered-R1_DNA1.Cyto_86-Scene-X15Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_60-Scene-X11Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_28-Scene-X11Y2_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_227-Scene-X9Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_50-Scene-X7Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_30-Scene-X8Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_226-Scene-X6Y9_c1c2.png', 'Registered-R1_DNA1.Cyto_138-Scene-X3Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_7-Scene-X5Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_191-Scene-X15Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_163-Scene-X1Y9_c1c2.png', 'Registered-R1_DNA1.Cyto_276-Scene-X9Y8_c1c2.png','Registered-R1_DNA1.Cyto_281-Scene-X8Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_256-Scene-X9Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_235-Scene-X12Y2_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_26-Scene-X16Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_129-Scene-X7Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_150-Scene-X5Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_190-Scene-X9Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_273-Scene-X6Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_175-Scene-X15Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_124-Scene-X10Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_223-Scene-X9Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_189-Scene-X5Y5_c1c2.png','Registered-R1_DNA1.Cyto_89-Scene-X4Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_52-Scene-X9Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_206-Scene-X2Y2_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_140-Scene-X10Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_184-Scene-X6Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_251-Scene-X16Y2_c1c2.png','Registered-R1_DNA1.Cyto_182-Scene-X1Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_180-Scene-X5Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_268-Scene-X6Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_178-Scene-X13Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_203-Scene-X8Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_72-Scene-X11Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_222-Scene-X8Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_75-Scene-X2Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_9-Scene-X6Y9_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_100-Scene-X15Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_251-Scene-X15Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_247-Scene-X13Y6_c1c2.png','Registered-R1_DNA1.Cyto_41-Scene-X11Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_282-Scene-X1Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_263-Scene-X2Y4_c1c2.png','Registered-R1_DNA1.Cyto_35-Scene-X16Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_115-Scene-X4Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_58-Scene-X16Y1_c1c2.png','Registered-R1_DNA1.Cyto_218-Scene-X9Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_139-Scene-X1Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_159-Scene-X12Y4_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_123-Scene-X8Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_79-Scene-X8Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_28-Scene-X12Y2_c1c2.png','Registered-R1_DNA1.Cyto_81-Scene-X12Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_166-Scene-X15Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_85-Scene-X10Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_241-Scene-X7Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_181-Scene-X6Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_58-Scene-X15Y1_c1c2.png','Registered-R1_DNA1.Cyto_272-Scene-X12Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_217-Scene-X1Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_1-Scene-X1Y8_c1c2.png',\n",
    " #'Registered-R1_DNA1.Cyto_6-Scene-X11Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_118-Scene-X13Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_193-Scene-X9Y7_c1c2.png','Registered-R1_DNA1.Cyto_188-Scene-X10Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_240-Scene-X10Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_272-Scene-X11Y3_c1c2.png','Registered-R1_DNA1.Cyto_94-Scene-X12Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_194-Scene-X11Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_146-Scene-X16Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_132-Scene-X13Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_114-Scene-X13Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_248-Scene-X8Y8_c1c2.png',\n",
    " 'Registered-R1_DNA1.Cyto_109-Scene-X10Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_45-Scene-X4Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_70-Scene-X5Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_202-Scene-X9Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_126-Scene-X14Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_74-Scene-X16Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_198-Scene-X7Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_153-Scene-X7Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_208-Scene-X15Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_154-Scene-X7Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_27-Scene-X6Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_36-Scene-X6Y1_c1c2.png',\n",
    " 'Registered-R1_DNA1.Cyto_48-Scene-X14Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_177-Scene-X16Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_157-Scene-X4Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_179-Scene-X13Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_200-Scene-X9Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_71-Scene-X8Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_147-Scene-X6Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_275-Scene-X6Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_161-Scene-X13Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_204-Scene-X16Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_253-Scene-X13Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_102-Scene-X1Y4_c1c2.png',\n",
    " 'Registered-R1_DNA1.Cyto_101-Scene-X10Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_24-Scene-X6Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_232-Scene-X7Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_131-Scene-X14Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_265-Scene-X14Y1_c1c2.png', 'Registered-R1_DNA1.Cyto_266-Scene-X7Y8_c1c2.png','Registered-R1_DNA1.Cyto_31-Scene-X16Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_234-Scene-X10Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_51-Scene-X1Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_103-Scene-X5Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_114-Scene-X14Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_63-Scene-X4Y9_c1c2.png',\n",
    " 'Registered-R1_DNA1.Cyto_73-Scene-X12Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_249-Scene-X3Y9_c1c2.png', 'Registered-R1_DNA1.Cyto_145-Scene-X1Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_125-Scene-X5Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_102-Scene-X2Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_120-Scene-X5Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_121-Scene-X11Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_87-Scene-X15Y4_c1c2.png', 'Registered-R1_DNA1.Cyto_110-Scene-X14Y8_c1c2.png', 'Registered-R1_DNA1.Cyto_212-Scene-X1Y9_c1c2.png', 'Registered-R1_DNA1.Cyto_25-Scene-X3Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_65-Scene-X14Y6_c1c2.png',\n",
    " 'Registered-R1_DNA1.Cyto_97-Scene-X16Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_221-Scene-X2Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_217-Scene-X2Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_52-Scene-X10Y2_c1c2.png', 'Registered-R1_DNA1.Cyto_149-Scene-X10Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_24-Scene-X5Y2_c1c2.png','Registered-R1_DNA1.Cyto_243-Scene-X13Y3_c1c2.png', 'Registered-R1_DNA1.Cyto_219-Scene-X8Y7_c1c2.png', 'Registered-R1_DNA1.Cyto_216-Scene-X1Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_76-Scene-X13Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_66-Scene-X15Y6_c1c2.png', 'Registered-R1_DNA1.Cyto_12-Scene-X12Y8_c1c2.png',\n",
    " 'Registered-R1_DNA1.Cyto_18-Scene-X13Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_93-Scene-X15Y5_c1c2.png', 'Registered-R1_DNA1.Cyto_240-Scene-X9Y2_c1c2.png','Registered-R1_DNA1.Cyto_187-Scene-X4Y2_c1c2.png'\n",
    " ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a56ab65",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#loop\n",
    "#labels = io.imread(f'{s_scene} nuclei{nuc_diam} - Nuclei Segmentation Basins.tif')\n",
    "#cell_labels = io.imread(f'{s_scene} cell{cell_diam} - Cell Segmentation Basins.tif')\n",
    "image_mpp = 1\n",
    "app = Mesmer()\n",
    "for s_file in ls_file:\n",
    "    img = io.imread(f'{segdir}/Mesmer_Images/{s_file}')\n",
    "    s_scene = s_file.split('_')[2]\n",
    "    print(s_scene)\n",
    "    a_test = np.expand_dims(img[:,:,0:2], axis=0)\n",
    "    rgb_test = create_rgb_image(a_test, channel_colors=['blue', 'green'])\n",
    "    #nuc\n",
    "    test_predictions = app.predict(a_test, image_mpp=image_mpp, compartment='nuclear')\n",
    "    labels = np.reshape(test_predictions[0],test_predictions[0].shape[0:2])\n",
    "    io.imsave(f'{segdir}/Mesmer_Segmentation/{s_scene} nuclei{image_mpp} - Nuclei Segmentation Basins.tif',labels,check_contrast=False)\n",
    "    #cell\n",
    "    cell_predictions = app.predict(a_test, image_mpp=image_mpp)\n",
    "    cell_labels = np.reshape(cell_predictions[0],cell_predictions[0].shape[0:2])\n",
    "    #clean up and save cells\n",
    "    df = pd.DataFrame(measure.regionprops_table(label_image=cell_labels,intensity_image=labels,properties=('label','area','mean_intensity')))\n",
    "    ls_drop = df[df.mean_intensity==0].label\n",
    "    cell_labels[np.isin(cell_labels, ls_drop)] = 0\n",
    "    io.imsave(f'{segdir}/Mesmer_Segmentation/{s_scene} cell{image_mpp} - Cell Segmentation Basins.tif',cell_labels,check_contrast=False)\n",
    "    #break"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12432317",
   "metadata": {},
   "source": [
    "## match"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8009565c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numba import jit\n",
    "import time\n",
    "from skimage import io, measure, segmentation\n",
    "#from deepcell.applications import Mesmer\n",
    "#from deepcell.utils.plot_utils import create_rgb_image\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "#oldcontainer,d_replace = segment.nuc_to_cell(labels,cell_labels)\n",
    "\n",
    "@jit(nopython=True, parallel=True)\n",
    "def _relabel_numba(ai_cell_labels, ai_cellid_relabel, o_dtype=np.uint32):\n",
    "    '''\n",
    "    version: 2021-12-00\n",
    "\n",
    "    input:\n",
    "        ai_cell_labels: (n,) one dimensional numpy array with cell labels.\n",
    "        ai_cellid_relabel: (n,2) dimensional numpy array which maps original cell_id to nuc_id.\n",
    "        o_dtype: data type inside the image array. standard is np.uint32,\n",
    "            because 32bit tiffs can be opened with most tiff compatible software.\n",
    "\n",
    "    output:\n",
    "        ai_cell_relabeld: (n,) one dimensional numpy array with cell labels that match the nucleus.\n",
    "\n",
    "    description:\n",
    "        use numba to quickly iterate over each label and replace pixels with new pixel values.\n",
    "    '''\n",
    "    # handle input\n",
    "    ai_cell_labels = ai_cell_labels.astype(o_dtype)\n",
    "    ai_cellid_relabel = ai_cellid_relabel.astype(o_dtype)\n",
    "\n",
    "    # do match\n",
    "    ai_cell_relabeld = ai_cell_labels.copy()\n",
    "    for i_cellid_cell, i_cellid_nuc in ai_cellid_relabel:\n",
    "        ai_cell_relabeld[ai_cell_labels == i_cellid_cell] = i_cellid_nuc\n",
    "\n",
    "    # output\n",
    "    return(ai_cell_relabeld)\n",
    "\n",
    "def _cell_to_nuc(ai_nuc_labels, ai_cell_labels):\n",
    "    '''\n",
    "    version: 2021-12-00\n",
    "\n",
    "    input:\n",
    "        ai_nuc_labels: (n,) one dimensional numpy array with nucleus labels.\n",
    "        ai_cell_labels: (n,) one dimensional numpy array with cell labels.\n",
    "\n",
    "    output:\n",
    "        ai_cellid_relabel: (n,2) dimensional numpy array which maps original cell_id to nuc_id.\n",
    "\n",
    "    description:\n",
    "        associate the largest nucleus contained in each cell segmentation.\n",
    "    '''\n",
    "    # get cell dominant nuclei label relationship\n",
    "    lli_cellid_relabel = []  # numba can't handle dictionary\n",
    "\n",
    "    # iterate over each cell label\n",
    "    for i_cellid in np.unique(ai_cell_labels):\n",
    "        if i_cellid == 0:\n",
    "            continue\n",
    "\n",
    "        # find all non-zero nuclei label contained within that cell mask\n",
    "        ai_nuc = ai_nuc_labels[ai_cell_labels == i_cellid]\n",
    "        ai_nuc = ai_nuc[ai_nuc != 0]\n",
    "        ai_nuc_unique = np.unique(ai_nuc)\n",
    "        # multiple nuclei, choose largest (most common pixels, i.e. mode)\n",
    "        if len(ai_nuc_unique) > 1:\n",
    "            ai_value, ai_count = np.unique(ai_nuc, return_counts=True) # does not work with numba\n",
    "            i_nucid_mode = ai_value[np.argmax(ai_count)]\n",
    "            if (i_cellid != i_nucid_mode):\n",
    "                lli_cellid_relabel.append([i_cellid, i_nucid_mode])\n",
    "        # one nucleus\n",
    "        elif len(ai_nuc_unique) == 1:\n",
    "            i_nucid_one = ai_nuc[0]\n",
    "            if (i_cellid != i_nucid_one):\n",
    "                lli_cellid_relabel.append([i_cellid, i_nucid_one])\n",
    "        # zero nuclei\n",
    "        else:\n",
    "            lli_cellid_relabel.append([i_cellid, 0])\n",
    "\n",
    "    # output\n",
    "    ai_cellid_relabel = np.array(lli_cellid_relabel) #, dtype=np.int64)\n",
    "    return(ai_cellid_relabel)\n",
    "\n",
    "def relabel_cells(labels,cell_labels):\n",
    "    start = time.time()\n",
    "    # ravel\n",
    "    ti_nuc_shape = labels.shape\n",
    "    ai_nuc_labels = labels.ravel()\n",
    "    ti_cell_shape = cell_labels.shape\n",
    "    ai_cell_labels = cell_labels.ravel()\n",
    "\n",
    "    # match cell id with nuclei id\n",
    "    print('get _cell_to_nuc_ listing ...')\n",
    "    ai_cellid_relabel = _cell_to_nuc(ai_nuc_labels=ai_nuc_labels, ai_cell_labels=ai_cell_labels)\n",
    "    \n",
    "    # set minimalistic variable type\n",
    "    print('running numba to relabel cell labels ...')\n",
    "    if (ai_cellid_relabel.max() < 2**8):\n",
    "        o_dtype = np.uint8\n",
    "    elif (ai_cellid_relabel.max() < 2**16):\n",
    "        o_dtype = np.uint16\n",
    "    elif (ai_cellid_relabel.max() < 2**32):\n",
    "        o_dtype = np.uint32\n",
    "    else:\n",
    "        # bue 20210810: currently, this will never happen. cellpose hardcoded the mask dtype to uint32.\n",
    "        # this limits cell indexing to: 2**32 -1 = 4'294'967'295 [cell].\n",
    "        o_dtype = np.uint64\n",
    "    # numba re-label cell id\n",
    "    ai_cell_relabeled = _relabel_numba(ai_cell_labels=ai_cell_labels, ai_cellid_relabel=ai_cellid_relabel, o_dtype=o_dtype)\n",
    "\n",
    "    # unravel\n",
    "    ai_nuc_labels.shape = ti_nuc_shape\n",
    "    ai_cell_relabeled.shape = ti_cell_shape\n",
    "    end = time.time()\n",
    "    print(f'time elapsed: {end - start}[sec]')\n",
    "    return(ai_cell_relabeled)\n",
    "\n",
    "#features\n",
    "def extract_feat(labels,intensity_image, properties=('centroid','mean_intensity','area','eccentricity')):\n",
    "    ''' \n",
    "    given labels and intensity image, extract features to dataframe\n",
    "    '''\n",
    "    props = measure.regionprops_table(labels,intensity_image, properties=properties)\n",
    "    df_prop = pd.DataFrame(props)\n",
    "    return(df_prop)\n",
    "\n",
    "def label_difference(labels,cell_labels):\n",
    "    '''\n",
    "    given matched nuclear and cell label IDs,return cell_labels minus labels\n",
    "    '''\n",
    "    overlap = cell_labels==labels\n",
    "    ring_rep = cell_labels.copy()\n",
    "    ring_rep[overlap] = 0\n",
    "    return(ring_rep)\n",
    "\n",
    "def index_format(df,s_scene,s_col='label'):\n",
    "    df.index = [f'{s_scene.split(\"-Scene-\")[0]}_scene{s_scene.split(\"-Scene-\")[1].split(\"_\")[0]}_cell{item}' for item in df.loc[:,s_col]]\n",
    "    return(df)\n",
    "\n",
    "def expand_label(labels,distance=3):\n",
    "    '''\n",
    "    expand the nucelar labels by a fixed number of pixels\n",
    "    '''\n",
    "    boundaries = segmentation.find_boundaries(labels,mode='outer') #thick\n",
    "    shrunk_labels = labels.copy()\n",
    "    shrunk_labels[boundaries] = 0\n",
    "    background = shrunk_labels == 0\n",
    "    distances, (i, j) = scipy.ndimage.distance_transform_edt(\n",
    "                background, return_indices=True\n",
    "            )\n",
    "\n",
    "    grown_labels = labels.copy()\n",
    "    mask = background & (distances <= distance)\n",
    "    grown_labels[mask] = shrunk_labels[i[mask], j[mask]]\n",
    "    ring_labels = grown_labels - shrunk_labels\n",
    "\n",
    "    return(ring_labels, grown_labels) #shrunk_labels, grown_labels,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56316537",
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_scene_done = ['10-Scene-X1Y5', '104-Scene-X8Y4', '107-Scene-X4Y7', '11-Scene-X13Y7', '113-Scene-X2Y5', '117-Scene-X13Y3',\n",
    " '117-Scene-X14Y3', '13-Scene-X5Y8', '131-Scene-X13Y4', '132-Scene-X14Y1', '134-Scene-X3Y8', '136-Scene-X9Y4', '14-Scene-X1Y8',\n",
    " '141-Scene-X11Y2', '143-Scene-X1Y2', '15-Scene-X1Y9', '16-Scene-X3Y4', '160-Scene-X13Y5', '162-Scene-X16Y7', '164-Scene-X9Y8',\n",
    " '17-Scene-X12Y4', '183-Scene-X3Y6', '195-Scene-X11Y6', '201-Scene-X14Y5', '213-Scene-X7Y7', '215-Scene-X4Y1', '22-Scene-X12Y6',\n",
    " '220-Scene-X10Y5', '229-Scene-X1Y7', '230-Scene-X11Y5', '233-Scene-X13Y7', '237-Scene-X16Y6', '239-Scene-X11Y3', '239-Scene-X12Y3',\n",
    " '241-Scene-X8Y3', '243-Scene-X14Y3', '246-Scene-X2Y6', '250-Scene-X2Y3', '252-Scene-X6Y7', '255-Scene-X12Y8', '259-Scene-X10Y3',\n",
    " '262-Scene-X10Y8', '263-Scene-X1Y4', '271-Scene-X11Y8', '280-Scene-X9Y5', '282-Scene-X2Y3', '284-Scene-X7Y2', '285-Scene-X7Y7',\n",
    " '3-Scene-X13Y6', '31-Scene-X15Y2', '34-Scene-X12Y5', '39-Scene-X1Y7', '4-Scene-X16Y3', '43-Scene-X5Y5',\n",
    " '46-Scene-X10Y7', '56-Scene-X2Y6', '57-Scene-X8Y6', '6-Scene-X12Y1', '68-Scene-X9Y1', '69-Scene-X5Y3', '70-Scene-X6Y3',\n",
    " '80-Scene-X4Y3', '83-Scene-X13Y2', '83-Scene-X14Y2', '95-Scene-X14Y2', '96-Scene-X9Y3']\n",
    "#from mplex_image import segment, features\n",
    "ls_scene = sorted(set([item.split(' ')[0] for item in os.listdir(f'{segdir}/Mesmer_Segmentation')])) #- set(ls_scene_done))\n",
    "len(ls_scene)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49060bec",
   "metadata": {},
   "outputs": [],
   "source": [
    "#matching #done\n",
    "'''\n",
    "for s_scene in ls_scene:\n",
    "    print(s_scene)\n",
    "    #s_scene = '46-Scene-X10Y7'\n",
    "    labels = io.imread(f'{segdir}/Mesmer_Segmentation/{s_scene} nuclei{image_mpp} - Nuclei Segmentation Basins.tif')\n",
    "    cell_labels = io.imread(f'{segdir}/Mesmer_Segmentation/{s_scene} cell{image_mpp} - Cell Segmentation Basins.tif')\n",
    "    cell_relabeled = relabel_cells(labels,cell_labels)\n",
    "    io.imsave(f'{segdir}/Mesmer_Segmentation/{s_scene} matchedcell{image_mpp} - Cell Segmentation Basins.tif',cell_relabeled)\n",
    "    #break\n",
    "print(labels.max())\n",
    "print(cell_labels.max())\n",
    "print(cell_relabeled.max())\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d1c3633",
   "metadata": {},
   "source": [
    "## extract features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a62a5a79",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(subdir)\n",
    "df_img = mpimage.parse_org()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c4d342e",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "image_mpp=1\n",
    "df_sample = pd.DataFrame()\n",
    "for s_scene in ls_scene:\n",
    "    labels = io.imread(f'{segdir}/Mesmer_Segmentation/{s_scene} nuclei{image_mpp} - Nuclei Segmentation Basins.tif')\n",
    "    cell_labels = io.imread(f'{segdir}/Mesmer_Segmentation/{s_scene} matchedcell{image_mpp} - Cell Segmentation Basins.tif')\n",
    "    s_dapi = f'Registered-R1_DNA1.Cyto_{s_scene}_c1c2.png'\n",
    "    dapi = io.imread(f'{segdir}/Mesmer_Images/{s_dapi}')[:,:,0]\n",
    "    #\n",
    "    df_feat = extract_feat(labels,dapi, properties=(['label','area','eccentricity','centroid']))\n",
    "    df_feat.columns = [f'{item}_nuclei' for item in df_feat.columns]\n",
    "    df_feat = index_format(df_feat,s_scene,'label_nuclei')   \n",
    "    cyto = label_difference(labels,cell_labels)\n",
    "    d_loc = {'nuclei':labels,'cell':cell_labels,'cytoplasm':cyto}\n",
    "    #markers (not duplicated\n",
    "    df_scene = df_img[df_img.scene==s_scene].sort_values('round_ord')\n",
    "    df_scene = df_scene[(df_scene.marker!='RutheniumTetroxide') & (~df_scene.index.str.contains('Registered-R11_HistoneH3'))& (~df_scene.index.str.contains('Registered-R24_HistoneH3'))]\n",
    "    #loop through markers\n",
    "    for s_index in df_scene.index:\n",
    "        df_marker = pd.DataFrame(index=df_feat.index)\n",
    "        s_marker = df_scene.loc[s_index,'marker']\n",
    "        print(f'extracting features {s_marker}')\n",
    "        if s_marker == 'DAPI':\n",
    "            s_marker = s_marker + f'{df_scene.loc[s_index,\"rounds\"].split(\"R\")[1]}'\n",
    "        intensity_image = io.imread(s_index)\n",
    "        for s_loc, a_loc in d_loc.items():\n",
    "            df_marker_loc = extract_feat(a_loc,intensity_image, properties=(['mean_intensity','label']))\n",
    "            df_marker_loc.columns = [f'{s_marker}_{s_loc}',f'label']\n",
    "            df_marker_loc = index_format(df_marker_loc,s_scene) ##  df_marker_loc.index]\n",
    "            df_marker_loc.drop(f'label',axis=1,inplace=True)\n",
    "            df_marker = df_marker.merge(df_marker_loc, left_index=True,right_index=True,how='left',suffixes=('',f'{s_marker}_{s_loc}'))\n",
    "        #extract nas\n",
    "        df_na = df_marker[df_marker.isna().any(axis=1)]\n",
    "        ls_na = [int(item.split('_cell')[1]) for item in df_na.index]\n",
    "        ring_labels, grown_labels = expand_label(labels,distance=3)\n",
    "        ring_labels[~np.isin(ring_labels,ls_na)] = 0\n",
    "        grown_labels[~np.isin(grown_labels,ls_na)] = 0\n",
    "        df_marker_ring = extract_feat(ring_labels,intensity_image, properties=(['mean_intensity','label']))\n",
    "        df_marker_grown = extract_feat(grown_labels,intensity_image, properties=(['mean_intensity','label']))\n",
    "        df_feat_na = df_marker_ring.merge(df_marker_grown,on='label',suffixes=('_ring','_grown'))\n",
    "        df_feat_na = index_format(df_feat_na,s_scene)\n",
    "        #fill in nas\n",
    "        #break\n",
    "        df_marker.loc[df_feat_na.index,f'{s_marker}_cytoplasm'] = df_feat_na.loc[:,f'mean_intensity_ring']\n",
    "        df_marker.loc[df_feat_na.index,f'{s_marker}_cell'] = df_feat_na.loc[:,f'mean_intensity_grown']   \n",
    "        df_feat = df_feat.merge(df_marker, left_index=True,right_index=True,how='left',suffixes=('',f'{s_marker}'))\n",
    "        #break\n",
    "    df_sample = df_sample.append(df_feat)\n",
    "    #break\n",
    "#df_sample.to_csv('features_IMC_MeanIntensity_Centroid_Shape.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0131362",
   "metadata": {},
   "source": [
    "## filter features "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42872161",
   "metadata": {},
   "outputs": [],
   "source": [
    "pwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1cfc6d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(segdir)\n",
    "df_mi = pd.read_csv(('features_IMC_MeanIntensity_Centroid_Shape.csv'),index_col=0)\n",
    "#df_mi = df_mi.drop(['DNA2_cell', 'DNA2_nuclei', 'DNA2_cytoplasm'],axis=1).dropna()\n",
    "#df_mi = df_mi.dropna() #618961 to 618811\n",
    "#150 cells are na in the markers: probaly too small for the expansion (get borders erased them)\n",
    "#df_mi.isna().sum().sort_values()[-3::].index\n",
    "#df_mi[df_mi.drop(['DNA2_cell', 'DNA2_nuclei', 'DNA2_cytoplasm'],axis=1).isna().ArgonDimers_cytoplasm].sort_values(by='area_nuclei')\n",
    "#df_mi.drop(['DNA2_cell', 'DNA2_nuclei', 'DNA2_cytoplasm'],axis=1).dropna()\n",
    "#imc segmentation had 855668/ 844498 w/o liver cells\n",
    "# mesmer seg has 618961 cells, 618811 NA free"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bff0e79",
   "metadata": {},
   "outputs": [],
   "source": [
    "#drop na df_norm is the watershed \n",
    "df_norm = pd.read_csv(f'{codedir}/data/20211221_BaselTMA_MeanIntensity.csv',index_col=0,low_memory=False)\n",
    "df_norm['slide_scene'] = [item.split('_')[2] + '_' + item.split('_')[3] for item in  df_norm.index]\n",
    "df_mi['slide_scene'] = [item.split('_cell')[0].replace('scene','') for item in  df_mi.index]\n",
    "df_norm = df_norm[~df_norm.slide_scene.isin(['Liver_X1Y1', 'Liver_X2Y1'])]\n",
    "print(set(df_norm.slide_scene) - set(df_mi.slide_scene))\n",
    "print(set(df_mi.slide_scene) - set(df_norm.slide_scene))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f89ed56",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "fig, ax = plt.subplots(figsize=(4,5),dpi=200)\n",
    "df_norm.groupby('slide_scene').Area.count().sort_values().plot(kind='barh',ax=ax,width=1)\n",
    "ax.yaxis.set_visible(False)\n",
    "ax.set_xlim(0,7000)\n",
    "ax.set_title('Watershed Segmentation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "404744f4",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(4,5),dpi=200)\n",
    "ls_index = df_norm.groupby('slide_scene').Area.count().sort_values().index\n",
    "df_mi.groupby('slide_scene').area_nuclei.count().loc[ls_index].plot(kind='barh',ax=ax,width=1)\n",
    "ax.yaxis.set_visible(False)\n",
    "ax.set_xlim(0,7000)\n",
    "ax.set_title('Mesmer Segmentation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d10a46a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_ws = df_norm.groupby('slide_scene').Area.count().sort_values()\n",
    "ls_index = df_norm.groupby('slide_scene').Area.count().sort_values().index\n",
    "df_mes = df_mi.groupby('slide_scene').area_nuclei.count().loc[ls_index]\n",
    "df_both = pd.DataFrame(data=df_ws).rename({'Area':'WS'},axis=1)\n",
    "df_both['MES'] = df_mes\n",
    "df_both.loc[ls_seg,'off'] = 'orange'\n",
    "df_both.fillna('blue',inplace=True)\n",
    "df_both.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd2fa968",
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "from scipy import stats\n",
    "fig, ax = plt.subplots(figsize=(4,3),dpi=300)\n",
    "sns.regplot(x=\"WS\", y=\"MES\", data=df_both,ax=ax,scatter_kws={'s':5,'alpha':0.5,\n",
    "                                                             'color':list(df_both['off'].values),'label':['__','b'],\n",
    "                                                            })\n",
    "r, pval= stats.pearsonr(df_both['WS'], df_both['MES'])\n",
    "ax.set_xlabel('Watershed cell count')\n",
    "ax.set_ylabel('Mesmer cell count')\n",
    "ax.legend(['__','outliers'],labelcolor='orange')\n",
    "ax.set_title(f'Watershed vs Mesmer\\nPearson r={r:.2} p={pval:.2}')\n",
    "print(pval)\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'{qcdir}/scatterplot_Mesmer_vs_Wastershed_cell_count.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "482ca9df",
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "df_both.loc[ls_seg_two,'two'] = 'orange'\n",
    "df_both.fillna('blue',inplace=True)\n",
    "fig, ax = plt.subplots(figsize=(4,3),dpi=300)\n",
    "sns.regplot(x=\"WS\", y=\"MES\", data=df_both,ax=ax,\n",
    "            scatter_kws={'s':5,'alpha':0.5,'color':list(df_both['two'].values),'label':['__','o']})\n",
    "r, pval= stats.pearsonr(df_both.loc[df_both.two=='blue','WS'], df_both.loc[df_both.two=='blue','MES'])\n",
    "ax.set_xlabel('Watershed cell count')\n",
    "ax.legend(['__','two WS ROIS'])\n",
    "ax.set_ylabel('Mesmer cell count')\n",
    "ax.set_title(f'Watershed vs Mesmer\\nPearson r={r:.2} p={pval:.2}')\n",
    "print(pval)\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'{qcdir}/scatterplot_Mesmer_vs_Wastershed_cell_count_two.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc8a7342",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(((df_ws - df_mes)/df_ws).mean())\n",
    "print(((df_ws - df_mes)/df_ws).median())\n",
    "print(len(df_ws))\n",
    "#16 have 50% less cells\n",
    "ls_seg = ((df_ws - df_mes)/df_ws).sort_values()[-16::].index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a40b2fdb",
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_seg = set(ls_seg) - set(ls_seg_two)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cabb298c",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(codedir)\n",
    "os.chdir('ome/masks')\n",
    "df_ws = pd.DataFrame(index=sorted(os.listdir()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "934d7abd",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(codedir)\n",
    "os.chdir('Segmentation/Mesmer_Segmentation')\n",
    "df_mes_all = pd.DataFrame(index=sorted(os.listdir()))\n",
    "df_mes = df_mes_all[df_mes_all.index.str.contains('matched')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1f557af",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_seg = pd.DataFrame(index=os.listdir(segdir))\n",
    "df_seg.drop(['Mesmer_Segmentation','features_IMC_MeanIntensity_Centroid_Shape.csv','test'],inplace=True)\n",
    "df_seg['slide'] = [item.split('-Scene-')[0].split('_')[-1] for item in df_seg.index]\n",
    "df_seg['scene'] =[item.split('-Scene-')[1].split('_')[0] for item in df_seg.index]\n",
    "df_seg['slide_scene'] = df_seg.slide + '_' + df_seg.scene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0081da0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mes['slide'] = [item.split('-Scene-')[0] for item in df_mes.index]\n",
    "df_mes['scene'] =[item.split('-Scene-')[1].split(' ')[0] for item in df_mes.index]\n",
    "df_mes['slide_scene'] = df_mes.slide + '_' + df_mes.scene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d115b549",
   "metadata": {},
   "outputs": [],
   "source": [
    "#BaselTMA_SP43_10.3kx13.7ky_6500x8000_6_20170928_29_19_X15Y2_216_a0_full_maks.tiff\n",
    "df_ws['slide'] = [item.split('_')[-6] for item in df_ws.index]\n",
    "df_ws['scene'] =[item.split('_')[-5] for item in df_ws.index]\n",
    "df_ws['slide_scene'] = df_ws.slide + '_' + df_ws.scene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a746f302",
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_seg_two = []\n",
    "for s_seg in sorted(set(df_seg.slide_scene)):\n",
    "    i_len = len(df_ws[df_ws.slide_scene==s_seg].index)\n",
    "    if i_len > 1:\n",
    "        ls_seg_two.append(s_seg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9afaac79",
   "metadata": {},
   "outputs": [],
   "source": [
    "for s_seg in ls_seg:\n",
    "    #s_seg = '104_X8Y4'\n",
    "    s_file_mes = df_mes[df_mes.slide_scene==s_seg].index[0]\n",
    "    #print(s_seg)\n",
    "    i_len = len(df_ws[df_ws.slide_scene==s_seg].index)\n",
    "    s_file_ws = df_ws[df_ws.slide_scene==s_seg].index[0]\n",
    "    s_file_img = df_seg[df_seg.slide_scene==s_seg].index[0]\n",
    "    matched = io.imread(f'{segdir}/Mesmer_Segmentation/{s_file_mes}')\n",
    "    wshed = io.imread(f'{codedir}/ome/masks/{s_file_ws}')\n",
    "    img = io.imread(f'{segdir}/{s_file_img}')\n",
    "    #img = \n",
    "    fig, ax = plt.subplots(1,3,figsize=(8,3),dpi=200)\n",
    "    img1 =io.imshow(matched,ax=ax[1],cmap='viridis')\n",
    "    #fig.colorbar(img1, ax=ax[0])\n",
    "    img2 = io.imshow(wshed,ax=ax[0],cmap='viridis')\n",
    "    #fig.colorbar(img2, ax=ax[1])\n",
    "    io.imshow(img,ax=ax[2])\n",
    "    ax[0].set_title(f'watershed: {wshed.max()}')\n",
    "    ax[1].set_title(f'Mesmer: {matched.max()}')\n",
    "    ax[2].set_title('Image')\n",
    "    ax[0].set_xlim(0,img.shape[1])\n",
    "    ax[0].set_ylim(img.shape[0],0)\n",
    "    ax[1].set_xlim(0,img.shape[1])\n",
    "    ax[1].set_ylim(img.shape[0],0)\n",
    "    \n",
    "    fig.suptitle(s_seg,y=0.9)\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'{qcdir}/segmentation_basins_{s_seg}.png')\n",
    "    if i_len > 1:\n",
    "        s_file_ws = df_ws[df_ws.slide_scene==s_seg].index[1]\n",
    "        wshed = io.imread(f'{codedir}/ome/masks/{s_file_ws}')\n",
    "        fig, ax = plt.subplots(figsize=(3,3),dpi=200)\n",
    "        io.imshow(wshed,ax=ax,cmap='viridis')\n",
    "    #break"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e6b923b",
   "metadata": {},
   "source": [
    "### testing seg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "309ae52b",
   "metadata": {},
   "outputs": [],
   "source": [
    "    #test_predictions = app.predict(a_test, image_mpp=image_mpp, compartment='nuclear')\n",
    "    #labels = np.reshape(test_predictions[0],test_predictions[0].shape[0:2])\n",
    "    #io.imsave(f'{segdir}/Mesmer_Segmentation/{s_scene} nuclei{image_mpp} - Nuclei Segmentation Basins.tif',labels,check_contrast=False)\n",
    "    #cell\n",
    "    #cell_predictions = app.predict(a_test, image_mpp=image_mpp)\n",
    "    #cell_labels = np.reshape(cell_predictions[0],cell_predictions[0].shape[0:2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1e7ae60",
   "metadata": {},
   "outputs": [],
   "source": [
    "io.imshow(labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eacd3100",
   "metadata": {},
   "outputs": [],
   "source": [
    "io.imshow(cell_labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76e812b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "#cleaned up\n",
    "io.imshow(cell_labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c273f5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(measure.regionprops_table(label_image=cell_labels,intensity_image=labels,properties=('label','area','mean_intensity')))\n",
    "ls_drop = df[df.mean_intensity==0].label\n",
    "cell_labels[np.isin(cell_labels, ls_drop)] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e3684d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#a_test = np.expand_dims(img[:,:,0:2], axis=0)\n",
    "#a_test.shape\n",
    "#rgb_test = create_rgb_image(a_test, channel_colors=['blue', 'green'])\n",
    "#whole cell\n",
    "#app = Mesmer()\n",
    "#test_predictions = app.predict(a_test, image_mpp=1)\n",
    "#io.imshow(test_predictions[0])\n",
    "#np.reshape(test_predictions[0],test_predictions[0].shape[0:2])\n",
    "#io.imsave(f'{segdir}/test/predictions.tiff',np.reshape(test_predictions[0],test_predictions[0].shape[0:2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5fe9b669",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "'''\n",
    "overlay_test = make_outline_overlay(rgb_data=rgb_test, predictions=test_predictions)\n",
    "idx=0\n",
    "i_size = 200\n",
    "fig, ax = plt.subplots(1, 2, figsize=(15, 15))\n",
    "ax[0].imshow(rgb_test[idx, 0:i_size,0:i_size,:])\n",
    "ax[1].imshow(overlay_test[idx, 0:i_size,0:i_size,:])\n",
    "ax[0].set_title('Raw data')\n",
    "ax[1].set_title('Predictions')\n",
    "plt.show()\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d5f8794",
   "metadata": {},
   "outputs": [],
   "source": [
    "#nuclear\n",
    "#test_predictions_nuc = app.predict(a_test, image_mpp=1, compartment='nuclear')\n",
    "#io.imshow(test_predictions_nuc[0])\n",
    "#np.reshape(test_predictions[0],test_predictions[0].shape[0:2])\n",
    "#io.imsave(f'{segdir}/test/predictions_nuc.tiff',np.reshape(test_predictions_nuc[0],test_predictions_nuc[0].shape[0:2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb0c9142",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "overlay_test_nuc = make_outline_overlay(rgb_data=rgb_test, predictions=test_predictions_nuc)\n",
    "idx=0\n",
    "i_size = 200\n",
    "fig, ax = plt.subplots(1, 2, figsize=(15, 15))\n",
    "ax[0].imshow(rgb_test[idx, 0:i_size,0:i_size,:])\n",
    "ax[1].imshow(overlay_test_nuc[idx, 0:i_size,0:i_size,:])\n",
    "ax[0].set_title('Raw data')\n",
    "ax[1].set_title('Predictions')\n",
    "plt.show()\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "081bbbb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "#napari clean up\n",
    "'''\n",
    "from skimage import measure\n",
    "df = pd.DataFrame(measure.regionprops_table(label_image=cell_labels,intensity_image=nuc_labels,properties=('label','area','mean_intensity')))\n",
    "ls_drop = df[df.mean_intensity==0].label\n",
    "cell_labels[np.isin(cell_labels, ls_drop)] = 0\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "974a5ea3",
   "metadata": {},
   "source": [
    "# image processing\n",
    "\n",
    "test differerent de-noising"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bc58f59",
   "metadata": {},
   "outputs": [],
   "source": [
    "# also tried simple\n",
    "# scipy.ndi.gaussian_filter()."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2da01adf",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(f'{segdir}/test')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9ec381c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#a_scale = skimage.exposure.rescale_intensity(dapi,in_range=(np.quantile(dapi,0.03),1.5*np.quantile(dapi,0.9999)))\n",
    "a_scale = dapi_scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00e52ed7",
   "metadata": {},
   "outputs": [],
   "source": [
    "io.imshow(a_scale[200:300,300:400])\n",
    "plt.title('Clip .9999 quantile')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e97e712b",
   "metadata": {},
   "outputs": [],
   "source": [
    "io.imsave('scaled.png',a_scale[200:600,300:700])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93603e06",
   "metadata": {},
   "outputs": [],
   "source": [
    "a_gamma = skimage.exposure.adjust_gamma(a_scale, gamma=0.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32664b29",
   "metadata": {},
   "outputs": [],
   "source": [
    "io.imshow(a_gamma[200:300,300:400])\n",
    "plt.title('Adjust gamma')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b4dcb7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "io.imsave('gamma.png',a_gamma[200:600,300:700])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b31b39db",
   "metadata": {},
   "outputs": [],
   "source": [
    "a_gauss = scipy.ndimage.gaussian_filter(a_gamma,sigma=1)\n",
    "io.imshow(a_gauss[200:300,300:400])\n",
    "plt.title('Gaussian blurr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81dee8b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "image = a_gamma[200:600,300:700]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cfd7539",
   "metadata": {},
   "outputs": [],
   "source": [
    "result_1 = unsharp_mask(image, radius=1, amount=1)\n",
    "result_2 = unsharp_mask(image, radius=5, amount=2)\n",
    "result_3 = unsharp_mask(image, radius=20, amount=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4039409",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=2, ncols=2,\n",
    "                         sharex=True, sharey=True, figsize=(10, 10))\n",
    "ax = axes.ravel()\n",
    "\n",
    "ax[0].imshow(image, cmap=plt.cm.gray)\n",
    "ax[0].set_title('Original image')\n",
    "ax[1].imshow(result_1, cmap=plt.cm.gray)\n",
    "ax[1].set_title('Enhanced image, radius=1, amount=1.0')\n",
    "ax[2].imshow(result_2, cmap=plt.cm.gray)\n",
    "ax[2].set_title('Enhanced image, radius=5, amount=2.0')\n",
    "ax[3].imshow(result_3, cmap=plt.cm.gray)\n",
    "ax[3].set_title('Enhanced image, radius=20, amount=1.0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b702be93",
   "metadata": {},
   "outputs": [],
   "source": [
    "io.imsave('result_1.png',(result_1*255).astype('uint8'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9c47ebd",
   "metadata": {},
   "outputs": [],
   "source": [
    "io.imsave('result_2.png',(result_2*255).astype('uint8'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb0473ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "io.imsave('result_3.png',(result_3*255).astype('uint8'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "add7ec79",
   "metadata": {},
   "outputs": [],
   "source": [
    "mip = np.stack(img_st).max(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "408cf667",
   "metadata": {},
   "outputs": [],
   "source": [
    "#added = np.clip((img_st[0] + img_st[1]),a_min=0,a_max=65535)\n",
    "#added_gamma = skimage.exposure.adjust_gamma(added, gamma=0.5)\n",
    "#pretty good\n",
    "noisy = a_gamma[200:600,300:700]\n",
    "tv = denoise_tv_chambolle(noisy, weight=0.25)\n",
    "#io.imsave('added_more_tv.png',(tv*255).astype('uint8'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0f7888b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#io.imshow(mip[200:600,300:700])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5dab2664",
   "metadata": {},
   "outputs": [],
   "source": [
    "#io.imsave('mip.png',mip[200:600,300:700])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ff698e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "#mip_gamma = skimage.exposure.adjust_gamma(mip, gamma=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7ebbcf1",
   "metadata": {},
   "outputs": [],
   "source": [
    "#io.imshow(mip_gamma[200:600,300:700])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d0cd1c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#io.imsave('gmip.png',mip_gamma[200:600,300:700])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "308550da",
   "metadata": {},
   "outputs": [],
   "source": [
    "#denoise\n",
    "noisy = a_gamma[200:300,300:400]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cb7fc99",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(8, 6), dpi=300,\n",
    "                       sharex=True, sharey=True)\n",
    "ax[0, 0].imshow(noisy,cmap='gray')\n",
    "ax[0, 0].axis('off')\n",
    "ax[0, 0].set_title('Noisy')\n",
    "ax[0, 1].imshow(denoise_tv_chambolle(noisy, weight=0.1),cmap='gray')\n",
    "ax[0, 1].axis('off')\n",
    "ax[0, 1].set_title('TV (0.1)')\n",
    "ax[0, 2].imshow(denoise_bilateral(noisy, sigma_color=0.05, sigma_spatial=15, \n",
    "                ),cmap='gray')\n",
    "ax[0, 2].axis('off')\n",
    "ax[0, 2].set_title('Bilateral (0.05)')\n",
    "ax[1, 0].imshow(denoise_wavelet(noisy,  rescale_sigma=True),cmap='gray')\n",
    "ax[1, 0].axis('off')\n",
    "ax[1, 0].set_title('Wavelet denoising')\n",
    "\n",
    "ax[1, 1].imshow(denoise_tv_chambolle(noisy, weight=0.2),cmap='gray')\n",
    "ax[1, 1].axis('off')\n",
    "ax[1, 1].set_title('TV (0.25)')\n",
    "ax[1, 2].imshow(denoise_bilateral(noisy, sigma_color=.2,sigma_spatial=15, \n",
    "                ),cmap='gray')\n",
    "ax[1, 2].axis('off')\n",
    "ax[1, 2].set_title('Bilateral (0.2)')\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "baa8012f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d945d5e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "noisy = added[200:250,300:350]\n",
    "fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(8, 6), dpi=300,\n",
    "                       sharex=True, sharey=True)\n",
    "ax[0, 0].imshow(noisy)\n",
    "ax[0, 0].axis('off')\n",
    "ax[0, 0].set_title('Noisy')\n",
    "ax[0, 1].imshow(denoise_tv_chambolle(noisy, weight=0.1))\n",
    "ax[0, 1].axis('off')\n",
    "ax[0, 1].set_title('TV')\n",
    "ax[0, 2].imshow(denoise_bilateral(noisy, sigma_color=0.05, sigma_spatial=15,\n",
    "                ))\n",
    "ax[0, 2].axis('off')\n",
    "ax[0, 2].set_title('Bilateral')\n",
    "ax[1, 0].imshow(denoise_wavelet(noisy,  rescale_sigma=True))\n",
    "ax[1, 0].axis('off')\n",
    "ax[1, 0].set_title('Wavelet denoising')\n",
    "\n",
    "ax[1, 1].imshow(denoise_tv_chambolle(noisy, weight=0.2))\n",
    "ax[1, 1].axis('off')\n",
    "ax[1, 1].set_title('(more) TV')\n",
    "ax[1, 2].imshow(denoise_bilateral(noisy, sigma_color=0.1, sigma_spatial=15,\n",
    "                ))\n",
    "ax[1, 2].axis('off')\n",
    "ax[1, 2].set_title('(more) Bilateral')\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4cbe765c",
   "metadata": {},
   "outputs": [],
   "source": [
    "added_gamma = skimage.exposure.adjust_gamma(added, gamma=0.5)\n",
    "#pretty good\n",
    "noisy = added_gamma[200:600,300:700]\n",
    "tv = denoise_tv_chambolle(noisy, weight=0.25)\n",
    "io.imsave('added_more_tv.png',(tv*255).astype('uint8'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c09198a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "noisy = added_gamma[200:250,300:350]\n",
    "fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(8, 6), dpi=300,\n",
    "                       sharex=True, sharey=True)\n",
    "ax[0, 0].imshow(noisy)\n",
    "ax[0, 0].axis('off')\n",
    "ax[0, 0].set_title('Noisy')\n",
    "ax[0, 1].imshow(denoise_tv_chambolle(noisy, weight=0.1))\n",
    "ax[0, 1].axis('off')\n",
    "ax[0, 1].set_title('TV')\n",
    "ax[0, 2].imshow(denoise_bilateral(noisy, sigma_color=0.05, sigma_spatial=15,\n",
    "                ))\n",
    "ax[0, 2].axis('off')\n",
    "ax[0, 2].set_title('Bilateral')\n",
    "ax[1, 0].imshow(denoise_wavelet(noisy,  rescale_sigma=True))\n",
    "ax[1, 0].axis('off')\n",
    "ax[1, 0].set_title('Wavelet denoising')\n",
    "\n",
    "ax[1, 1].imshow(denoise_tv_chambolle(noisy, weight=0.25))\n",
    "ax[1, 1].axis('off')\n",
    "ax[1, 1].set_title('(more) TV')\n",
    "ax[1, 2].imshow(denoise_bilateral(noisy, sigma_color=0.1, sigma_spatial=15,\n",
    "                ))\n",
    "ax[1, 2].axis('off')\n",
    "ax[1, 2].set_title('(more) Bilateral')\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ed347e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#pretty good\n",
    "noisy = added_gamma[200:600,300:700]\n",
    "tv = denoise_tv_chambolle(noisy, weight=0.25)\n",
    "io.imsave('added_more_tv.png',(tv*255).astype('uint8'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85584959",
   "metadata": {},
   "outputs": [],
   "source": [
    "io.imshow(noisy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43ded57e",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "footprint = morphology.disk(2)\n",
    "res = morphology.white_tophat(noisy, footprint)\n",
    "fig, ax = plt.subplots(ncols=3, figsize=(10, 5),dpi=300)\n",
    "ax[0].set_title('Original')\n",
    "ax[0].imshow(noisy, cmap='gray')\n",
    "ax[1].set_title('White tophat')\n",
    "ax[1].imshow(res, cmap='gray')\n",
    "ax[2].set_title('Complementary')\n",
    "ax[2].imshow(noisy - res, cmap='gray')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e696e61",
   "metadata": {},
   "outputs": [],
   "source": [
    "tv = denoise_tv_chambolle(noisy, weight=0.25)\n",
    "footprint = morphology.disk(2)\n",
    "res = morphology.white_tophat(tv, footprint)\n",
    "fig, ax = plt.subplots(ncols=3, figsize=(10, 5),dpi=300)\n",
    "ax[0].set_title('Original')\n",
    "ax[0].imshow(tv, cmap='gray')\n",
    "ax[1].set_title('White tophat')\n",
    "ax[1].imshow(res, cmap='gray')\n",
    "ax[2].set_title('Complementary')\n",
    "ax[2].imshow(tv - res, cmap='gray')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c264fbff",
   "metadata": {},
   "outputs": [],
   "source": [
    "noisy = added_gamma[200:300,300:400]\n",
    "image = denoise_tv_chambolle(noisy, weight=0.25)\n",
    "result_1 = unsharp_mask(image, radius=1, amount=1)\n",
    "result_2 = unsharp_mask(image, radius=5, amount=1)\n",
    "result_3 = unsharp_mask(image, radius=20, amount=1)\n",
    "\n",
    "fig, axes = plt.subplots(nrows=2, ncols=2,\n",
    "                         sharex=True, sharey=True, figsize=(10, 10))\n",
    "ax = axes.ravel()\n",
    "\n",
    "ax[0].imshow(image, cmap=plt.cm.gray)\n",
    "ax[0].set_title('Original image')\n",
    "ax[1].imshow(result_1, cmap=plt.cm.gray)\n",
    "ax[1].set_title('Enhanced image, radius=1, amount=1.0')\n",
    "ax[2].imshow(result_2, cmap=plt.cm.gray)\n",
    "ax[2].set_title('Enhanced image, radius=5, amount=1.0')\n",
    "ax[3].imshow(result_3, cmap=plt.cm.gray)\n",
    "ax[3].set_title('Enhanced image, radius=20, amount=1.0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6f4a0de",
   "metadata": {},
   "outputs": [],
   "source": [
    "noisy = added_gamma[200:300,300:400]\n",
    "image = denoise_tv_chambolle(noisy, weight=0.25)\n",
    "result_1 = unsharp_mask(image, radius=1, amount=1)\n",
    "result_2 = unsharp_mask(image, radius=1, amount=2)\n",
    "result_3 = unsharp_mask(image, radius=1, amount=4)\n",
    "\n",
    "fig, axes = plt.subplots(nrows=2, ncols=2,\n",
    "                         sharex=True, sharey=True, figsize=(10, 10))\n",
    "ax = axes.ravel()\n",
    "\n",
    "ax[0].imshow(image, cmap=plt.cm.gray)\n",
    "ax[0].set_title('Original image')\n",
    "ax[1].imshow(result_1, cmap=plt.cm.gray)\n",
    "ax[1].set_title('Enhanced image, radius=1, amount=1.0')\n",
    "ax[2].imshow(result_2, cmap=plt.cm.gray)\n",
    "ax[2].set_title('Enhanced image, radius=1, amount=2.0')\n",
    "ax[3].imshow(result_3, cmap=plt.cm.gray)\n",
    "ax[3].set_title('Enhanced image, radius=1, amount=4.0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c743375",
   "metadata": {},
   "outputs": [],
   "source": [
    "#pretty good\n",
    "noisy = added_gamma[200:600,300:700]\n",
    "tv = denoise_tv_chambolle(noisy, weight=0.25)\n",
    "unsharp = unsharp_mask(tv, radius=1, amount=4)\n",
    "io.imsave('added_more_tv_unsharp.png',(unsharp*255).astype('uint8'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13cc39be",
   "metadata": {},
   "outputs": [],
   "source": [
    "#pretty bad\n",
    "\n",
    "tvu = denoise_tv_chambolle(unsharp, weight=0.1)\n",
    "io.imsave('added_more_tv_us_tv.png',(tvu*255).astype('uint8'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "899c9d63",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python3.9.5",
   "language": "python",
   "name": "python3.9.5"
  },
  "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}