{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using FISSA with CNMF from MATLAB\n",
    "\n",
    "CNMF is blind source separation toolbox for cell detection and signal extraction. \n",
    "\n",
    "Here we illustrate how one can use the ROIs detected by CNMF, and use FISSA to extract and decontaminate the traces.\n",
    "\n",
    "In this tutorial, we assume the user is using the [MATLAB implementation of CNMF](https://github.com/flatironinstitute/CaImAn-MATLAB). As such, this also serves as a tutorial on how to import data from MATLAB into Python to use with FISSA.\n",
    "\n",
    "However, note that there is also a [Python implementation of CNMF](https://github.com/flatironinstitute/CaImAn), which you can use instead to keep your whole workflow in Python.\n",
    "\n",
    "**Reference:**\n",
    "Pnevmatikakis, E.A., Soudry, D., Gao, Y., Machado, T., Merel, J., ... and Paninski, L. Simultaneous denoising, deconvolution, and demixing of calcium imaging data. *Neuron*, **89**(2):285-299, 2016. doi: [10.1016/j.neuron.2015.11.037](http://dx.doi.org/10.1016/j.neuron.2015.11.037)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FISSA package\n",
    "import fissa\n",
    "\n",
    "# For plotting our results, import numpy and matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Need this utility from scipy to load data from matfiles\n",
    "from scipy.io import loadmat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running CNMF in MATLAB, and importing into Python\n",
    "\n",
    "We ran CNMF in MATLAB using the [run_pipeline.m](https://github.com/flatironinstitute/CaImAn-MATLAB/blob/master/run_pipeline.m) script available from the CNMF repository on our example data (found at [../exampleData/20150529/](https://github.com/rochefort-lab/fissa/tree/1.0.0/examples/exampleData)).\n",
    "\n",
    "We [saved](https://mathworks.com/help/matlab/ref/save.html) the `Coor` and `F_df` variables generated by that script into a `.mat` file ([cNMFdata.mat](https://github.com/rochefort-lab/fissa/blob/1.0.0/examples/cNMFdata.mat?raw=true)) which we now load here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data from cNMFdata.mat file\n",
    "cNMFdata = loadmat(\"cNMFdata\")[\"dat\"]\n",
    "\n",
    "# Get the F_df recording traces out of the loaded object\n",
    "F_df = cNMFdata[\"F_df\"][0, 0]\n",
    "\n",
    "# Get the ROI outlines out of the loaded object\n",
    "Coor = cNMFdata[\"Coor\"][0, 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Show detected cells\n",
    "\n",
    "Let's render the ROIs using matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting lines surrounding each of the ROIs\n",
    "plt.figure(figsize=(7, 7))\n",
    "\n",
    "for i_cell in range(len(Coor)):\n",
    "    x = Coor[i_cell, 0][0]\n",
    "    y = Coor[i_cell, 0][1]\n",
    "    # Plot border around cells\n",
    "    plt.plot(x, y)\n",
    "\n",
    "# Invert the y-axis because image co-ordinates are labelled from top-left\n",
    "plt.gca().invert_yaxis()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running FISSA on cells detected by CNMF\n",
    "\n",
    "FISSA needs ROIs to be provided either as an ImageJ zip file, or a set of numpy arrays.\n",
    "\n",
    "CNMF can output ROIs in coordinates (as we imported above), which can be directly read into FISSA.\n",
    "A given ROI after importing from MATLAB is given as\n",
    "\n",
    "```python\n",
    "Coor[i, 0]\n",
    "```\n",
    "\n",
    "FISSA expects a set of rois to be given as a list of lists,\n",
    "```python\n",
    "[[roiA1, roiA2, roiA3, ...]]\n",
    "```\n",
    "so we will need to change the format of the ROIs first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numROI = len(Coor)\n",
    "rois_FISSA = [[Coor[i, 0][0], Coor[i, 0][1]] for i in range(numROI)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which can then be put into FISSA and run as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_folder = \"fissa_cnmf_example\"\n",
    "tiff_folder = \"exampleData/20150529/\"\n",
    "\n",
    "experiment = fissa.Experiment(tiff_folder, [rois_FISSA], output_folder)\n",
    "experiment.separate(redo_prep=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting the results\n",
    "\n",
    "Let's plot the traces for ROIs as they were detected by CNMF, and after removing neuropile with FISSA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fetch the colormap object for Cynthia Brewer's Paired color scheme\n",
    "cmap = plt.get_cmap(\"Paired\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select which trial (TIFF index) to plot\n",
    "trial = 0\n",
    "\n",
    "# Plot the mean image and ROIs from the FISSA experiment\n",
    "plt.figure(figsize=(7, 7))\n",
    "plt.imshow(experiment.means[trial], cmap=\"gray\")\n",
    "\n",
    "XLIM = plt.xlim()\n",
    "YLIM = plt.ylim()\n",
    "\n",
    "for i_roi in range(len(experiment.roi_polys)):\n",
    "    # Plot border around ROI\n",
    "    for contour in experiment.roi_polys[i_roi, trial][0]:\n",
    "        plt.plot(\n",
    "            contour[:, 1],\n",
    "            contour[:, 0],\n",
    "            color=cmap((i_roi * 2 + 1) % cmap.N),\n",
    "        )\n",
    "\n",
    "# ROI co-ordinates are half a pixel outside the image,\n",
    "# so we reset the axis limits\n",
    "plt.xlim(XLIM)\n",
    "plt.ylim(YLIM)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot all ROIs and trials\n",
    "\n",
    "# Get the number of ROIs and trials\n",
    "n_roi = experiment.result.shape[0]\n",
    "n_trial = experiment.result.shape[1]\n",
    "\n",
    "# Find the maximum signal intensities for each ROI\n",
    "roi_max_raw = [\n",
    "    np.max([np.max(experiment.raw[i_roi, i_trial][0]) for i_trial in range(n_trial)])\n",
    "    for i_roi in range(n_roi)\n",
    "]\n",
    "roi_max_result = [\n",
    "    np.max([np.max(experiment.result[i_roi, i_trial][0]) for i_trial in range(n_trial)])\n",
    "    for i_roi in range(n_roi)\n",
    "]\n",
    "roi_max = np.maximum(roi_max_raw, roi_max_result)\n",
    "\n",
    "# Plot our figure using subplot panels\n",
    "plt.figure(figsize=(16, 10))\n",
    "for i_roi in range(n_roi):\n",
    "    for i_trial in range(n_trial):\n",
    "        # Make subplot axes\n",
    "        i_subplot = 1 + i_trial * n_roi + i_roi\n",
    "        plt.subplot(n_trial, n_roi, i_subplot)\n",
    "        # Plot the data\n",
    "        plt.plot(\n",
    "            experiment.raw[i_roi][i_trial][0, :],\n",
    "            label=\"Raw (CNMF)\",\n",
    "            color=cmap((i_roi * 2) % cmap.N),\n",
    "        )\n",
    "        plt.plot(\n",
    "            experiment.result[i_roi][i_trial][0, :],\n",
    "            label=\"FISSA\",\n",
    "            color=cmap((i_roi * 2 + 1) % cmap.N),\n",
    "        )\n",
    "        # Labels and boiler plate\n",
    "        plt.ylim([-0.05 * roi_max[i_roi], roi_max[i_roi] * 1.05])\n",
    "        if i_roi == 0:\n",
    "            plt.ylabel(\n",
    "                \"Trial {}\\n\\nSignal intensity\\n(candela per unit area)\".format(\n",
    "                    i_trial + 1\n",
    "                )\n",
    "            )\n",
    "        if i_trial == 0:\n",
    "            plt.title(\"ROI {}\".format(i_roi))\n",
    "        if i_trial == n_trial - 1:\n",
    "            plt.xlabel(\"Time (frame number)\")\n",
    "            plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The figure shows the raw signal from the ROI identified by CNMF (pale), and after decontaminating with FISSA (dark).\n",
    "The hues match the ROI locations drawn above.\n",
    "Each column shows the results from one of the ROIs detected by CNMF.\n",
    "Each row shows the results from one of the three trials."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python",
   "language": "python",
   "name": "python"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}