napari_3d_image_processing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Napari demo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import packages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Click on **[this link](../desktop)** to open a new window where napari will be visible. Only required for mybinder\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ['DISPLAY'] = ':1.0'\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "import napari\n",
    "import skimage\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import requests\n",
    "import skimage.feature\n",
    "import pandas as pd\n",
    "import trackpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%gui qt5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Download the data from Zenodo."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This dataset is a time-lapse of a zebrafish embry acquired by confocal microscopy. See [here](https://zenodo.org/record/1211599#.XZSnAOf7TOR) for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.isfile('Sample/cxcr4aMO2_290112.lsm'):\n",
    "    os.makedirs('Sample/',exist_ok=True)\n",
    "    myfile = requests.get('https://zenodo.org/record/1211599/files/cxcr4aMO2_290112.lsm?download=1', allow_redirects=True)\n",
    "    open('Sample/cxcr4aMO2_290112.lsm', 'wb').write(myfile.content)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Image import and plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We rescale the image in z by a factor of 5 to clearly show all features. We only keep the later time points for the demo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "##USE THIS ON MYBINDER (RESIZED IMAGE AND FEWER TIME POINTS\n",
    "'''image_import = skimage.external.tifffile.TiffFile('Sample/cxcr4aMO2_290112.lsm') \n",
    "time_min = 100\n",
    "time_max = 130\n",
    "image_time = np.zeros((time_max-time_min,22,256,256))\n",
    "for t in range(time_min, time_max):\n",
    "    image_time[t-time_min,:,:,:] = image_import.asarray(key=2*np.arange(t*44,(t+1)*44))[::2,::2,::2]\n",
    "'''  \n",
    "\n",
    "##USE THIS ON A REGULAR COMPUTER\n",
    "image_time = skimage.io.imread('Sample/cxcr4aMO2_290112.lsm', \n",
    "                               plugin='tifffile')\n",
    "image_time = image_time[0,100:150,:,:,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "view3 = napari.Viewer(ndisplay=3)\n",
    "view3.add_image(np.swapaxes(image_time[:,:,:,:],1,2),scale=[1,5,1],\n",
    "                colormap ='green',blending='additive');\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nuclei detection\n",
    "We use a simple gaussian filtering followed by local maxima detection to find all the nuclei and create an array with time, x, y,z positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_dots = []\n",
    "filtered = np.zeros(image_time.shape)\n",
    "for t in range(image_time.shape[0]):\n",
    "    imG = skimage.filters.gaussian(image_time[t,:,:,:],sigma =2,\n",
    "                               preserve_range=True)\n",
    "    filtered[t,:,:,:] = imG\n",
    "    locmax = skimage.feature.peak_local_max(imG,min_distance=2,\n",
    "                                            threshold_abs=20)\n",
    "    newloc = np.stack((t*np.ones(locmax.shape[0]),\n",
    "                       locmax[:,1],5*locmax[:,0],locmax[:,2]),axis =1)\n",
    "    all_dots.append(newloc)\n",
    "spots = np.concatenate(all_dots) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The array containing spot detections can simply be plotted in Napari. If the dimensions are in the right order, napari naturally plots them correctly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "view3 = napari.Viewer(ndisplay=3)\n",
    "view3.add_image(np.swapaxes(filtered[:,:,:,:],1,2),scale=[1,5,1],\n",
    "                colormap ='cyan',blending='additive');\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Points layer 'Points' at 0x13b9637f0>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "view3.add_points(spots,face_color='blue',opacity=0.5,size=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tracking\n",
    "\n",
    "Finally we use [trackpy](https://soft-matter.github.io/trackpy/v0.3.2/) to track the different nuclei. For that we just have to turn our time,x,y,z array into a Pandas dataframe and then call trackpy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Frame 49: 219 trajectories present.\n"
     ]
    }
   ],
   "source": [
    "spots_pd = pd.DataFrame(spots,columns=['frame','x','y','z'])\n",
    "traj = trackpy.link_df(spots_pd,search_range=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We just exclude very short tracks to make the plot more readable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_traj = []\n",
    "for ind in traj.particle.unique():\n",
    "    curr_traj = traj[traj.particle==ind][['x','y','z']].values\n",
    "    if curr_traj.shape[0]>5:\n",
    "        all_traj.append(curr_traj)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally we add the trajectories to the plot as lines:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Shapes layer 'Shapes' at 0x13fce66d8>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "view3.add_shapes(all_traj, shape_type='path',edge_color='red')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}