Gad1 ALM Sample 1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sample ID: HCR 3.11 S5  \n",
    "Section: 200um coronal section from ALM  \n",
    "Mouse: ANM471122 C57BI/6NCrl, sex M, DOB 10/21/2019  \n",
    "Genes labelled: Gad1 (514 channel), Ndnf (546 channel), Vip (594 channel), Sst (647 channel)  \n",
    "Microscope: Zeiss 880 (regular confocal mode)  \n",
    "Objective: 40X oil immersion "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-02T06:01:30.013111Z",
     "start_time": "2020-04-02T06:01:29.998111Z"
    },
    "init_cell": true,
    "run_control": {
     "marked": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline    \n",
    "\n",
    "import napari\n",
    "from PIL import Image\n",
    "\n",
    "import numpy as np\n",
    "import cv2\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import time\n",
    "\n",
    "from os.path import sep\n",
    "import pickle as pkl\n",
    "import importlib\n",
    "\n",
    "from utils import diff_gauss\n",
    "from utils import find_threshold\n",
    "\n",
    "from segmentation import get_masks\n",
    "from segmentation import get_background\n",
    "\n",
    "from plotting import *\n",
    "\n",
    "from lipofuscin import annotated_lipo_rois\n",
    "from lipofuscin import filter_lipo\n",
    "\n",
    "from expression import *\n",
    "\n",
    "from classification import *\n",
    "\n",
    "from cell_depths import *\n",
    "\n",
    "# Enable autoreload of local modules\n",
    "%load_ext autoreload\n",
    "%autoreload 2        # Autoreload all modules except those excluded using %aimport\n",
    "%aimport napari, PIL.Image, numpy, cv2, matplotlib.pyplot, time, importlib \n",
    "\n",
    "# from mpl_toolkits import mplot3d\n",
    "# import matplotlib.colors as colors\n",
    "# import matplotlib.path as mpltpath\n",
    "\n",
    "\n",
    "# import h5py\n",
    "\n",
    "# from bisect import bisect_left\n",
    "\n",
    "# from skimage.measure import profile_line\n",
    "# from skimage.filters import gaussian\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Specify data path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-02T00:36:07.549671Z",
     "start_time": "2020-04-02T00:36:07.544671Z"
    },
    "hidden": true,
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "data_path = 'G:\\\\Shared drives\\\\as_share\\\\HCR\\\\HCR 3.11\\\\S5'\n",
    "metadata_file = 'metadata.pkl'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-02T05:18:19.264306Z",
     "start_time": "2020-04-02T05:18:19.210306Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of planes: 10\n"
     ]
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "\n",
    "um_per_px = 0.13\n",
    "raw_image_path = '{0}{1}S5_raw_tiff_images'.format(data_path, sep)\n",
    "plane_nos = range(1, 11)\n",
    "base_filename = 'S5__z'\n",
    "genes = ['Gad1', 'Ndnf', 'Sst', 'Vip']\n",
    "channel_names = {'Gad1': '_ChS1-T3_ORG',\n",
    "                 'Ndnf': '_ChS2-T1_ORG',\n",
    "                 'Vip': '_ChS2-T3_ORG',\n",
    "                 'Sst': '_ChS2-T2_ORG'}\n",
    "n_planes = len(plane_nos)\n",
    "print('Number of planes: {0}'.format(n_planes))\n",
    "\n",
    "# Create 4D array to store images\n",
    "img = Image.open('{0}{1}{2}{3}{4}.tif'.format(raw_image_path, sep, base_filename, str(plane_nos[0]).zfill(2),\n",
    "                                             channel_names[genes[0]]))\n",
    "\n",
    "h = img.height\n",
    "w = img.width\n",
    "\n",
    "cell_data_file = 'cell_data.pkl'\n",
    "cell_pixels_file = 'cell_pixels.pkl'\n",
    "all_cell_pixels_file = 'all_cell_pixels.pkl'\n",
    "bg_pixels_file = 'bg_pixels.pkl'\n",
    "all_bg_pixels_file = 'all_bg_pixels.pkl'\n",
    "lipo_rois_file = 'lipo_rois_cells.pkl'\n",
    "lipo_pixels_cells_file = 'lipo_pixels_cells.pkl'\n",
    "\n",
    "cell_filt_values_file = 'cell_filt_values.pkl'\n",
    "bg_filt_values_file = 'bg_filt_values.pkl'\n",
    "anno_lipo_filt_values_file = 'anno_lipo_filt_values.pkl'\n",
    "\n",
    "seg_image_path = '{0}{1}S5_blurred_tiff_images'.format(data_path, sep)\n",
    "plots_path = '{0}{1}Plots'.format(data_path, sep)\n",
    "filt_image_path = '{0}{1}Filtered_tiff_images'.format(data_path, sep)\n",
    "\n",
    "metadata = {}\n",
    "metadata['raw_image_path'] = raw_image_path\n",
    "metadata['filt_image_path'] = filt_image_path\n",
    "metadata['plane_nos'] = plane_nos\n",
    "metadata['base_filename'] = base_filename\n",
    "metadata['genes'] = genes\n",
    "metadata['channel_names'] = channel_names\n",
    "metadata['n_planes'] = n_planes\n",
    "metadata['h'] = h\n",
    "metadata['w'] = w\n",
    "metadata['um_per_px'] = um_per_px\n",
    "\n",
    "metadata['cell_data_file'] = cell_data_file\n",
    "metadata['cell_pixels_file'] = cell_pixels_file\n",
    "metadata['bg_pixels_file'] = bg_pixels_file\n",
    "metadata['all_bg_pixels_file'] = all_bg_pixels_file\n",
    "metadata['all_cell_pixels_file'] = all_cell_pixels_file\n",
    "metadata['lipo_rois_file'] = lipo_rois_file\n",
    "metadata['lipo_pixels_cells_file'] = lipo_pixels_cells_file\n",
    "\n",
    "metadata['cell_filt_values_file'] = cell_filt_values_file\n",
    "metadata['bg_filt_values_file'] = bg_filt_values_file\n",
    "metadata['anno_lipo_filt_values_file'] = anno_lipo_filt_values_file\n",
    "    \n",
    "metadata['seg_image_path'] = seg_image_path\n",
    "metadata['plots_path'] = plots_path\n",
    "    \n",
    "with open('{0}{1}{2}'.format(data_path, sep, metadata_file), 'wb') as f:\n",
    "    pkl.dump(metadata, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Segmentation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Get cell pixels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-02T04:28:22.425502Z",
     "start_time": "2020-04-02T03:08:00.501693Z"
    },
    "hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "357 cell masks found\n",
      "Finding pixels in cells\n",
      "Cell 0: 38.0 seconds\n",
      "Cell 2: 61.0 seconds\n",
      "Cell 4: 84.0 seconds\n",
      "Cell 6: 107.0 seconds\n",
      "Cell 8: 122.0 seconds\n",
      "Cell 10: 144.0 seconds\n",
      "Cell 12: 189.0 seconds\n",
      "Cell 14: 204.0 seconds\n",
      "Cell 16: 218.0 seconds\n",
      "Cell 18: 240.0 seconds\n",
      "Cell 20: 263.0 seconds\n",
      "Cell 22: 292.0 seconds\n",
      "Cell 24: 315.0 seconds\n",
      "Cell 26: 337.0 seconds\n",
      "Cell 28: 374.0 seconds\n",
      "Cell 30: 411.0 seconds\n",
      "Cell 32: 441.0 seconds\n",
      "Cell 34: 456.0 seconds\n",
      "Cell 36: 485.0 seconds\n",
      "Cell 38: 522.0 seconds\n",
      "Cell 40: 537.0 seconds\n",
      "Cell 42: 552.0 seconds\n",
      "Cell 44: 574.0 seconds\n",
      "Cell 46: 619.0 seconds\n",
      "Cell 48: 663.0 seconds\n",
      "Cell 50: 686.0 seconds\n",
      "Cell 52: 700.0 seconds\n",
      "Cell 54: 723.0 seconds\n",
      "Cell 56: 745.0 seconds\n",
      "Cell 58: 768.0 seconds\n",
      "Cell 60: 791.0 seconds\n",
      "Cell 62: 806.0 seconds\n",
      "Cell 64: 836.0 seconds\n",
      "Cell 66: 858.0 seconds\n",
      "Cell 68: 881.0 seconds\n",
      "Cell 70: 896.0 seconds\n",
      "Cell 72: 926.0 seconds\n",
      "Cell 74: 949.0 seconds\n",
      "Cell 76: 972.0 seconds\n",
      "Cell 78: 995.0 seconds\n",
      "Cell 80: 1026.0 seconds\n",
      "Cell 82: 1049.0 seconds\n",
      "Cell 84: 1071.0 seconds\n",
      "Cell 86: 1101.0 seconds\n",
      "Cell 88: 1123.0 seconds\n",
      "Cell 90: 1138.0 seconds\n",
      "Cell 92: 1168.0 seconds\n",
      "Cell 94: 1183.0 seconds\n",
      "Cell 96: 1220.0 seconds\n",
      "Cell 98: 1242.0 seconds\n",
      "Cell 100: 1272.0 seconds\n",
      "Cell 102: 1309.0 seconds\n",
      "Cell 104: 1331.0 seconds\n",
      "Cell 106: 1346.0 seconds\n",
      "Cell 108: 1368.0 seconds\n",
      "Cell 110: 1398.0 seconds\n",
      "Cell 112: 1420.0 seconds\n",
      "Cell 114: 1452.0 seconds\n",
      "Cell 116: 1491.0 seconds\n",
      "Cell 118: 1522.0 seconds\n",
      "Cell 120: 1556.0 seconds\n",
      "Cell 122: 1587.0 seconds\n",
      "Cell 124: 1618.0 seconds\n",
      "Cell 126: 1641.0 seconds\n",
      "Cell 128: 1672.0 seconds\n",
      "Cell 130: 1718.0 seconds\n",
      "Cell 132: 1749.0 seconds\n",
      "Cell 134: 1779.0 seconds\n",
      "Cell 136: 1810.0 seconds\n",
      "Cell 138: 1840.0 seconds\n",
      "Cell 140: 1878.0 seconds\n",
      "Cell 142: 1893.0 seconds\n",
      "Cell 144: 1923.0 seconds\n",
      "Cell 146: 1945.0 seconds\n",
      "Cell 148: 1968.0 seconds\n",
      "Cell 150: 1998.0 seconds\n",
      "Cell 152: 2021.0 seconds\n",
      "Cell 154: 2058.0 seconds\n",
      "Cell 156: 2081.0 seconds\n",
      "Cell 158: 2096.0 seconds\n",
      "Cell 160: 2126.0 seconds\n",
      "Cell 162: 2157.0 seconds\n",
      "Cell 164: 2194.0 seconds\n",
      "Cell 166: 2224.0 seconds\n",
      "Cell 168: 2271.0 seconds\n",
      "Cell 170: 2300.0 seconds\n",
      "Cell 172: 2330.0 seconds\n",
      "Cell 174: 2353.0 seconds\n",
      "Cell 176: 2383.0 seconds\n",
      "Cell 178: 2413.0 seconds\n",
      "Cell 180: 2443.0 seconds\n",
      "Cell 182: 2467.0 seconds\n",
      "Cell 184: 2497.0 seconds\n",
      "Cell 186: 2527.0 seconds\n",
      "Cell 188: 2557.0 seconds\n",
      "Cell 190: 2573.0 seconds\n",
      "Cell 192: 2603.0 seconds\n",
      "Cell 194: 2634.0 seconds\n",
      "Cell 196: 2658.0 seconds\n",
      "Cell 198: 2689.0 seconds\n",
      "Cell 200: 2719.0 seconds\n",
      "Cell 202: 2742.0 seconds\n",
      "Cell 204: 2780.0 seconds\n",
      "Cell 206: 2810.0 seconds\n",
      "Cell 208: 2841.0 seconds\n",
      "Cell 210: 2879.0 seconds\n",
      "Cell 212: 2917.0 seconds\n",
      "Cell 214: 2940.0 seconds\n",
      "Cell 216: 2971.0 seconds\n",
      "Cell 218: 3000.0 seconds\n",
      "Cell 220: 3023.0 seconds\n",
      "Cell 222: 3046.0 seconds\n",
      "Cell 224: 3069.0 seconds\n",
      "Cell 226: 3099.0 seconds\n",
      "Cell 228: 3129.0 seconds\n",
      "Cell 230: 3144.0 seconds\n",
      "Cell 232: 3174.0 seconds\n",
      "Cell 234: 3204.0 seconds\n",
      "Cell 236: 3227.0 seconds\n",
      "Cell 238: 3241.0 seconds\n",
      "Cell 240: 3257.0 seconds\n",
      "Cell 242: 3287.0 seconds\n",
      "Cell 244: 3324.0 seconds\n",
      "Cell 246: 3339.0 seconds\n",
      "Cell 248: 3354.0 seconds\n",
      "Cell 250: 3383.0 seconds\n",
      "Cell 252: 3405.0 seconds\n",
      "Cell 254: 3436.0 seconds\n",
      "Cell 256: 3450.0 seconds\n",
      "Cell 258: 3472.0 seconds\n",
      "Cell 260: 3502.0 seconds\n",
      "Cell 262: 3524.0 seconds\n",
      "Cell 264: 3555.0 seconds\n",
      "Cell 266: 3584.0 seconds\n",
      "Cell 268: 3607.0 seconds\n",
      "Cell 270: 3636.0 seconds\n",
      "Cell 272: 3651.0 seconds\n",
      "Cell 274: 3674.0 seconds\n",
      "Cell 276: 3711.0 seconds\n",
      "Cell 278: 3741.0 seconds\n",
      "Cell 280: 3778.0 seconds\n",
      "Cell 282: 3808.0 seconds\n",
      "Cell 284: 3830.0 seconds\n",
      "Cell 286: 3860.0 seconds\n",
      "Cell 288: 3882.0 seconds\n",
      "Cell 290: 3912.0 seconds\n",
      "Cell 292: 3941.0 seconds\n",
      "Cell 294: 3964.0 seconds\n",
      "Cell 296: 3993.0 seconds\n",
      "Cell 298: 4016.0 seconds\n",
      "Cell 300: 4038.0 seconds\n",
      "Cell 302: 4060.0 seconds\n",
      "Cell 304: 4083.0 seconds\n",
      "Cell 306: 4112.0 seconds\n",
      "Cell 308: 4135.0 seconds\n",
      "Cell 310: 4157.0 seconds\n",
      "Cell 312: 4202.0 seconds\n",
      "Cell 314: 4232.0 seconds\n",
      "Cell 316: 4270.0 seconds\n",
      "Cell 318: 4285.0 seconds\n",
      "Cell 320: 4299.0 seconds\n",
      "Cell 322: 4329.0 seconds\n",
      "Cell 324: 4352.0 seconds\n",
      "Cell 326: 4381.0 seconds\n",
      "Cell 328: 4411.0 seconds\n",
      "Cell 330: 4434.0 seconds\n",
      "Cell 332: 4463.0 seconds\n",
      "Cell 334: 4493.0 seconds\n",
      "Cell 336: 4523.0 seconds\n",
      "Cell 338: 4553.0 seconds\n",
      "Cell 340: 4575.0 seconds\n",
      "Cell 342: 4598.0 seconds\n",
      "Cell 344: 4635.0 seconds\n",
      "Cell 346: 4658.0 seconds\n",
      "Cell 348: 4687.0 seconds\n",
      "Cell 350: 4709.0 seconds\n",
      "Cell 352: 4755.0 seconds\n",
      "Cell 354: 4777.0 seconds\n",
      "Cell 356: 4807.0 seconds\n",
      "Cell pixels found in G:\\Shared drives\\as_share\\HCR\\HCR 3.11\\S5\\all_cell_pixels.pkl\n"
     ]
    }
   ],
   "source": [
    "cell_data, n_cells = get_masks.get_masks(data_path, metadata_file)\n",
    "cell_pixels = get_masks.get_cell_pixels(data_path, metadata_file)\n",
    "all_cell_pixels = get_masks.get_all_cell_pixels(data_path, metadata_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get background pixels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-02T06:15:39.353112Z",
     "start_time": "2020-04-02T06:15:25.806112Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finding pixels in local background of cells\n"
     ]
    }
   ],
   "source": [
    "bg_pixels = get_background.get_bg_pixels(data_path, metadata_file)\n",
    "all_bg_pixels = get_background.get_all_bg_pixels(data_path, metadata_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Segment cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lipofuscin exclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Annotated lipofuscin rois"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-02T03:07:02.821693Z",
     "start_time": "2020-04-02T03:07:02.806693Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14 rois\n",
      "Data loaded\n"
     ]
    }
   ],
   "source": [
    "l_rois = annotated_lipo_rois.get_rois(data_path, metadata_file)\n",
    "all_lipo_pixels, lipo_pixels_roi = annotated_lipo_rois.get_cell_lipo_pixels(data_path, metadata_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Detect lipofuscin rois in cells and background"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-02T11:01:15.489551Z",
     "start_time": "2020-04-02T10:59:34.992351Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Filtering Ndnf images: 0.0 seconds\n",
      "Plane 1, 0.0 seconds\n",
      "Saving plane 0: 11.0 seconds\n",
      "Plane 2, 11.0 seconds\n",
      "Saving plane 1: 31.0 seconds\n",
      "Plane 3, 31.0 seconds\n",
      "Saving plane 2: 58.0 seconds\n",
      "Plane 4, 58.0 seconds\n",
      "Saving plane 3: 84.0 seconds\n",
      "Plane 5, 85.0 seconds\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "Can not convert from object to float64.",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-126-1787c2ad3ca2>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      1\u001b[0m cell_filt_values, bg_filt_values, annotated_lipo_filt_values = filter_lipo.filter(data_path, metadata_file, \n\u001b[0;32m      2\u001b[0m                                                                                    \u001b[0mgenes\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;34m'Ndnf'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m                                                                                   save_tiffs = True)\n\u001b[0m",
      "\u001b[1;32mG:\\Shared drives\\as_share\\HCR\\FISH Image Analysis\\lipofuscin\\filter_lipo.py\u001b[0m in \u001b[0;36mfilter\u001b[1;34m(data_path, metadata_file, genes, sigma_small, sigma_large, thresh_scale, save_tiffs)\u001b[0m\n\u001b[0;32m     66\u001b[0m                                                             channel_names[gene]))\n\u001b[0;32m     67\u001b[0m             \u001b[0mim_array\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mimg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 68\u001b[1;33m             \u001b[0mim_diff\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdiff_gauss\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdiff_gauss\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msigma_small\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msigma_large\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mim_array\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     69\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     70\u001b[0m             \u001b[0mpx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mall_cell_pixels\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mall_cell_pixels\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m==\u001b[0m \u001b[0mplane\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mG:\\Shared drives\\as_share\\HCR\\FISH Image Analysis\\utils\\diff_gauss.py\u001b[0m in \u001b[0;36mdiff_gauss\u001b[1;34m(sigma_small, sigma_large, img, plot_kernel)\u001b[0m\n\u001b[0;32m     27\u001b[0m     \u001b[0mim_diff\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mimg\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     28\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 29\u001b[1;33m     \u001b[0mim_filt_large\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgaussian\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mimg\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msigma\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msigma_large\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     30\u001b[0m     \u001b[0mim_filt_small\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgaussian\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mimg\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msigma\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msigma_small\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     31\u001b[0m     \u001b[0mim_diff\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mim_filt_small\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mim_filt_large\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\AppData\\Local\\Continuum\\anaconda3\\envs\\fish\\lib\\site-packages\\skimage\\filters\\_gaussian.py\u001b[0m in \u001b[0;36mgaussian\u001b[1;34m(image, sigma, output, mode, cval, multichannel, preserve_range, truncate)\u001b[0m\n\u001b[0;32m    113\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msigma\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m!=\u001b[0m \u001b[0mimage\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mndim\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    114\u001b[0m             \u001b[0msigma\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mconcatenate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msigma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 115\u001b[1;33m     \u001b[0mimage\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mconvert_to_float\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mimage\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpreserve_range\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    116\u001b[0m     return ndi.gaussian_filter(image, sigma, mode=mode, cval=cval,\n\u001b[0;32m    117\u001b[0m                                truncate=truncate)\n",
      "\u001b[1;32m~\\AppData\\Local\\Continuum\\anaconda3\\envs\\fish\\lib\\site-packages\\skimage\\_shared\\utils.py\u001b[0m in \u001b[0;36mconvert_to_float\u001b[1;34m(image, preserve_range)\u001b[0m\n\u001b[0;32m    254\u001b[0m             \u001b[0mimage\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mimage\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfloat\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    255\u001b[0m     \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 256\u001b[1;33m         \u001b[0mimage\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mimg_as_float\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mimage\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    257\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mimage\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\AppData\\Local\\Continuum\\anaconda3\\envs\\fish\\lib\\site-packages\\skimage\\util\\dtype.py\u001b[0m in \u001b[0;36mimg_as_float\u001b[1;34m(image, force_copy)\u001b[0m\n\u001b[0;32m    428\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    429\u001b[0m     \"\"\"\n\u001b[1;32m--> 430\u001b[1;33m     \u001b[1;32mreturn\u001b[0m \u001b[0mconvert\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mimage\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfloating\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mforce_copy\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    431\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    432\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\AppData\\Local\\Continuum\\anaconda3\\envs\\fish\\lib\\site-packages\\skimage\\util\\dtype.py\u001b[0m in \u001b[0;36mconvert\u001b[1;34m(image, dtype, force_copy, uniform)\u001b[0m\n\u001b[0;32m    248\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mdtype_in\u001b[0m \u001b[1;32min\u001b[0m \u001b[0m_supported_types\u001b[0m \u001b[1;32mand\u001b[0m \u001b[0mdtype_out\u001b[0m \u001b[1;32min\u001b[0m \u001b[0m_supported_types\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    249\u001b[0m         raise ValueError(\"Can not convert from {} to {}.\"\n\u001b[1;32m--> 250\u001b[1;33m                          .format(dtypeobj_in, dtypeobj_out))\n\u001b[0m\u001b[0;32m    251\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    252\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0mkind_in\u001b[0m \u001b[1;32min\u001b[0m \u001b[1;34m'ui'\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mValueError\u001b[0m: Can not convert from object to float64."
     ]
    }
   ],
   "source": [
    "# Filter images \n",
    "cell_filt_values, bg_filt_values, annotated_lipo_filt_values = filter_lipo.filter(data_path, metadata_file, \n",
    "                                                                                   genes = ['Ndnf'],\n",
    "                                                                                  save_tiffs = True)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reload modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-02T10:59:24.186551Z",
     "start_time": "2020-04-02T10:59:23.995351Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<module 'utils.find_threshold' from 'G:\\\\Shared drives\\\\as_share\\\\HCR\\\\FISH Image Analysis\\\\utils\\\\find_threshold.py'>"
      ]
     },
     "execution_count": 125,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Ideally, modules should be auto-reloaded \n",
    "importlib.reload(get_masks)\n",
    "importlib.reload(get_background)\n",
    "\n",
    "importlib.reload(annotated_lipo_rois)\n",
    "importlib.reload(filter_lipo)\n",
    "\n",
    "importlib.reload(diff_gauss)\n",
    "importlib.reload(find_threshold)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Debugging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-02T10:52:56.398351Z",
     "start_time": "2020-04-02T10:52:45.823951Z"
    },
    "run_control": {
     "marked": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Filtering Gad1 images: 0.0 seconds\n",
      "Plane 1, 0.0 seconds\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'save_tiffs' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-110-78724d9e24d1>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     73\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     74\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 75\u001b[1;33m         \u001b[1;32mif\u001b[0m \u001b[0msave_tiffs\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     76\u001b[0m             \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Saving plane {0}: {1} seconds'\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mplane\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mround\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtime\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mt0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     77\u001b[0m             imwrite('{0}{1}{2}_plane{3}_filt_sigma={4},{5}_thresh={6}.tif'.format(filt_image_path, sep, gene,\n",
      "\u001b[1;31mNameError\u001b[0m: name 'save_tiffs' is not defined"
     ]
    }
   ],
   "source": [
    "sigma_small = 5\n",
    "sigma_large = 10\n",
    "thresh_scale = 1\n",
    "save_tiffs = True\n",
    "\n",
    "import time\n",
    "import pickle as pkl\n",
    "from os.path import sep\n",
    "from PIL import Image\n",
    "from utils import diff_gauss\n",
    "from utils import find_threshold\n",
    "from tifffile import imwrite\n",
    "from segmentation import get_masks\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "# Load metadata\n",
    "with open('{0}{1}{2}'.format(data_path, sep, metadata_file), 'rb') as f:\n",
    "    metadata = pkl.load(f)\n",
    "raw_image_path = metadata['raw_image_path']\n",
    "filt_image_path = metadata['filt_image_path']\n",
    "n_planes = metadata['n_planes']\n",
    "h = metadata['h']\n",
    "w = metadata['w']\n",
    "base_filename = metadata['base_filename']\n",
    "channel_names = metadata['channel_names']\n",
    "if genes == None:\n",
    "    genes = metadata['genes']\n",
    "\n",
    "# Load cell pixels\n",
    "all_cell_pixels_file = metadata['all_cell_pixels_file']\n",
    "with open('{0}{1}{2}'.format(data_path, sep, all_cell_pixels_file), 'rb') as f:\n",
    "    all_cell_pixels = pkl.load(f)\n",
    "\n",
    "cell_pixels_file = metadata['cell_pixels_file']\n",
    "with open('{0}{1}{2}'.format(data_path, sep, cell_pixels_file), 'rb') as f:\n",
    "    cell_pixels = pkl.load(f)\n",
    "\n",
    "# Load background pixels\n",
    "bg_pixels_file = metadata['bg_pixels_file']\n",
    "with open('{0}{1}{2}'.format(data_path, sep, bg_pixels_file), 'rb') as f:\n",
    "    bg_pixels = pkl.load(f)\n",
    "\n",
    "# Load annotated lipofuscin pixels\n",
    "lipo_pixels_cells_file = metadata['lipo_pixels_cells_file']\n",
    "with open('{0}{1}{2}'.format(data_path, sep, lipo_pixels_cells_file), 'rb') as f:\n",
    "    data_dict = pkl.load(f)\n",
    "all_lipo_pixels = data_dict['all_lipo_pixels']\n",
    "\n",
    "# Allocate variables to store pixel values in filtered images\n",
    "cell_filt_values = {cell: {} for cell in cell_pixels.keys()}\n",
    "bg_filt_values = {cell: {} for cell in cell_pixels.keys()}\n",
    "annotated_lipo_filt_values = []\n",
    "\n",
    "t0 = time.time()\n",
    "im_thresh = np.zeros([h, w])\n",
    "for gene in genes:\n",
    "\n",
    "    print('Filtering {0} images: {1} seconds'.format(gene, np.round(time.time() - t0)))\n",
    "\n",
    "    for plane in range(n_planes):\n",
    "        print('Plane {0}, {1} seconds'.format(plane + 1, np.round(time.time() - t0)))\n",
    "        img = Image.open('{0}{1}{2}{3}{4}.tif'.format(raw_image_path, sep, base_filename, str(plane + 1).zfill(2),\n",
    "                                                        channel_names[gene]))\n",
    "        im_array = np.array(img)\n",
    "        im_diff = diff_gauss.diff_gauss(sigma_small, sigma_large, im_array)\n",
    "\n",
    "        px = all_cell_pixels[all_cell_pixels[:, 0] == plane, 1:]\n",
    "        points = im_diff[px[:, 0], px[:, 1]]\n",
    "        thresh = find_threshold.find_threshold(points, thresh_scale)\n",
    "\n",
    "        ix = np.where(im_diff > thresh)\n",
    "        im_thresh[ix] = im_diff[ix]\n",
    "\n",
    "\n",
    "        if save_tiffs:\n",
    "            print('Saving plane {0}: {1} seconds'.format(plane, np.round(time.time() - t0)))\n",
    "            imwrite('{0}{1}{2}_plane{3}_filt_sigma={4},{5}_thresh={6}.tif'.format(filt_image_path, sep, gene,\n",
    "                                                                            str(plane + 1).zfill(2),\n",
    "                                                                            sigma_small, sigma_large,\n",
    "                                                                            np.round(thresh)),\n",
    "                                                                            im_thresh.astype(np.int16))\n",
    "\n",
    "        # Get pixel values for cells and background in filtered images\n",
    "        cells_in_plane = get_masks.get_cells_in_plane(data_path, metadata_file, plane)\n",
    "        for cell in cells_in_plane:\n",
    "\n",
    "            cell_px = cell_pixels[cell][plane]\n",
    "            cell_filt_values[cell][plane] = im_thresh[cell_px[0], cell_px[1]]\n",
    "\n",
    "            bg_px = bg_pixels[cell][plane]\n",
    "            bg_filt_values[cell][plane] = im_thresh[bg_px[:, 0], bg_px[:, 1]]\n",
    "\n",
    "        # Get pixel values for annotated lipofuscin pixels in filtered images\n",
    "        lipo_px = all_lipo_pixels[all_lipo_pixels[:, 0] == plane, 1:]\n",
    "        annotated_lipo_filt_values = np.append(annotated_lipo_filt_values, im_thresh[lipo_px[:, 0], lipo_px[:, 0]])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Initialization Cell",
  "kernelspec": {
   "display_name": "fish",
   "language": "python",
   "name": "fish"
  },
  "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.7.6"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "oldHeight": 122.66666599999999,
   "position": {
    "height": "635px",
    "left": "1366px",
    "right": "20px",
    "top": "122px",
    "width": "309px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "varInspector_section_display": "block",
   "window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}