exchange_aa_working.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns; sns.set()\n",
    "from matplotlib import pyplot as plt, cm, colors\n",
    "import matplotlib.gridspec as gridspec\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "sns.set(style=\"darkgrid\")\n",
    "sns.set(style=\"darkgrid\")\n",
    "from matplotlib import pyplot as plt, cm, colors\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "import glob\n",
    "import os.path as ospath\n",
    "import os\n",
    "from sys import executable\n",
    "from subprocess import check_output\n",
    "from PyQt5.QtWidgets import QFileDialog, QApplication\n",
    "from IPython.display import HTML\n",
    "\n",
    "from scipy import optimize\n",
    "from scipy.spatial import distance\n",
    "from scipy import linalg\n",
    "from scipy import signal\n",
    "from sklearn.cluster import MeanShift, estimate_bandwidth\n",
    "\n",
    "from picasso.picasso import io\n",
    "from picasso.picasso.postprocess import link, compute_dark_times\n",
    "from picasso.picasso.render import render\n",
    "from picasso.picasso.gui.render import estimate_kinetic_rate, fit_cum_exp\n",
    "\n",
    "# define visualization options \n",
    "%matplotlib inline\n",
    "%gui qt\n",
    "\n",
    "# define colors:\n",
    "blue = \"#4C72B0\"\n",
    "orange = \"#DD8452\"\n",
    "red = \"#C44E52\"\n",
    "gray = \"#90A8CE\"\n",
    "\n",
    "\n",
    "def OpenFileDialog():\n",
    "    file = check_output([executable, __file__])\n",
    "    return file.strip()\n",
    "\n",
    "\n",
    "def gui_fname(dir=None):\n",
    "    \"\"\"\n",
    "    Select a file via a dialog and return the file name.\n",
    "    \"\"\"\n",
    "    if dir is None: \n",
    "        dir =\"./\"\n",
    "\n",
    "    app = QApplication([dir])\n",
    "    fname = QFileDialog.getExistingDirectory(None, \"Select a folder...\", \n",
    "            dir)\n",
    "    if isinstance(fname, tuple):\n",
    "        return fname[0]\n",
    "    else: \n",
    "        return str(fname)\n",
    "\n",
    "\n",
    "def load_files(dirname):\n",
    "    \n",
    "    os.chdir(dirname)\n",
    "    files = glob.glob(\"*.hdf5\")\n",
    "    \n",
    "    if files:\n",
    "        print(\"{} HDF5 files found.\".format(len(files)))\n",
    "    else:\n",
    "        print(\"No HDF5 files found at: {}\".format(dirname))\n",
    "            \n",
    "    return files\n",
    "\n",
    "\n",
    "def load_data(path):\n",
    "\n",
    "    data, info = io.load_locs(path)\n",
    "    try:\n",
    "        pixelsize = info[1][\"Pixelsize\"]\n",
    "    except:  \n",
    "        print(\"No pixelsize found in yaml file. Default 130 nm used.\")\n",
    "  \n",
    "\n",
    "    # convert px to nm\n",
    "    data.x *= pixelsize\n",
    "    data.y *= pixelsize\n",
    "    data.x_pick_rot *= pixelsize\n",
    "    data.y_pick_rot *= pixelsize\n",
    "    \n",
    "    return data, info, pixelsize\n",
    "\n",
    "\n",
    "def double_gaus(x,a,x0,sigma, b, x1, sigma1):\n",
    "    return a*np.exp(-(x-x0)**2/(2*sigma**2)) + b*np.exp(-(x-x1)**2/(2*sigma1**2))\n",
    "\n",
    "def gaus(x,a,x0,sigma):\n",
    "    return a*np.exp(-(x-x0)**2/(2*sigma**2)) \n",
    "\n",
    "\n",
    "def find_peaks(data, binning=100, axes=\"y\"):\n",
    "    \n",
    "    if axes == \"y\":\n",
    "        column = \"y_pick_rot\"\n",
    "    elif axes == \"x\":\n",
    "        column = \"x_pick_rot\"\n",
    "    elif axes == \"xyz\":\n",
    "        column = 2\n",
    "    \n",
    "    # find peaks\n",
    "    bandwidth = estimate_bandwidth(data[column].reshape(-1, 1), quantile=0.2, n_samples=binning)\n",
    "    #print(\"estimated bandwidth: \"+str(bandwidth))\n",
    "    ms = MeanShift(bandwidth=bandwidth, bin_seeding=True)\n",
    "    ms.fit(data[column].reshape(-1, 1))\n",
    "    labels = ms.labels_\n",
    "    peaks = np.sort(ms.cluster_centers_[0:2],axis=None) # assuming that there are two large peaks from the two rings\n",
    "    peak1 = np.float(peaks[0])\n",
    "    peak2 = np.float(peaks[1])\n",
    "    estimated_peaks = {0:peak1, 1:peak2}\n",
    "\n",
    "    # histogram\n",
    "    n, bins = np.histogram(data[column], bins=binning)\n",
    "    bins = bins[1:]\n",
    "    hist_data = [n, bins]\n",
    "\n",
    "    # fit peaks\n",
    "    p0 = [peak1/2, peak1, 40, peak2/2, peak2, 40]\n",
    "    try:\n",
    "        p_fit, p_cov = optimize.curve_fit(double_gaus, bins, n, p0=p0)\n",
    "    except:\n",
    "        p_fit = [0,0,0,1000,1000,1000]\n",
    "        \n",
    "    \n",
    "    # check order of fitted peaks (peak 1 < peak 2)\n",
    "    if p_fit[1] > p_fit[4]:\n",
    "        p_temp = p_fit.copy()\n",
    "        p_fit[0:3]=p_temp[3:6]\n",
    "        p_fit[3:6]=p_temp[0:3]\n",
    "       \n",
    "        \n",
    "    # p_fit: amplitued_1, center_1, width_1, amplitude_2, center_2, width_2\n",
    "        \n",
    "    return estimated_peaks, p_fit, hist_data\n",
    "\n",
    "def find_peak(data, binning=100, axes=\"y\"):\n",
    "    \n",
    "    if axes == \"y\":\n",
    "        column = \"y_pick_rot\"\n",
    "    elif axes == \"x\":\n",
    "        column = \"x_pick_rot\"\n",
    "    elif axes == \"xyz\":\n",
    "        column = 2\n",
    "    \n",
    "    # find peak\n",
    "    bandwidth = estimate_bandwidth(data[column].reshape(-1, 1), quantile=0.2, n_samples=binning)\n",
    "    #print(\"estimated bandwidth: \"+str(bandwidth))\n",
    "    ms = MeanShift(bandwidth=bandwidth, bin_seeding=True)\n",
    "    ms.fit(data[column].reshape(-1, 1))\n",
    "    labels = ms.labels_\n",
    "   # print(ms.cluster_centers_[0:2])\n",
    "    peaks = np.sort(ms.cluster_centers_[0:2], axis=None)\n",
    "    peak1 = np.float(ms.cluster_centers_[0])\n",
    "    estimated_peak = {0:peak1}\n",
    "\n",
    "    # histogram\n",
    "    n, bins = np.histogram(data[column], bins=binning)\n",
    "    bins = bins[1:]\n",
    "    hist_data = [n, bins]\n",
    "\n",
    "    # fit peaks\n",
    "    p0 = [peak1/2, peak1, 40]\n",
    "    try:\n",
    "        p_fit, p_cov = optimize.curve_fit(gaus, bins, n, p0=p0)\n",
    "    except:\n",
    "        p_fit = [0,0,0]\n",
    "        \n",
    "    return estimated_peak, p_fit, hist_data\n",
    "\n",
    "\n",
    "def plot_peak_dist(data, hist_data, peaks, p_fit, binning=100,cutoff=1, axes=\"y\", ax=None):\n",
    "    \n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "        \n",
    "    if axes == \"y\":\n",
    "        column = \"y_pick_rot\"\n",
    "    elif axes == \"x\":\n",
    "        column = \"x_pick_rot\"\n",
    "    elif axes == \"xyz\":\n",
    "        column == 1\n",
    "        \n",
    "    peak1 = peaks[0]\n",
    "    peak2 = peaks[1]\n",
    "    \n",
    "    n = hist_data[0]\n",
    "    bins = hist_data[1]\n",
    "    \n",
    "    #fig, ax = plt.subplots(figsize=(9, 5))\n",
    "    binwidth = (max(bins)-min(bins)) / binning\n",
    "    ax.bar(bins, n, width=binwidth, color=gray)\n",
    "    xlin = np.linspace(0, data[column].max(), 1000)\n",
    "    ax.plot(xlin, gaus(xlin,*p_fit[0:3]), c=red, linewidth=2)\n",
    "    ax.plot(xlin, gaus(xlin,*p_fit[3:6]), c=red, linewidth=2)\n",
    "    ax.axvline(p_fit[1]-(p_fit[2]*cutoff),c=blue, linewidth=2, linestyle=\"--\")\n",
    "    ax.axvline(p_fit[1]+(p_fit[2]*cutoff),c=blue, linewidth=2, linestyle=\"--\")\n",
    "    ax.axvline(p_fit[4]-(p_fit[5]*cutoff),c=orange, linewidth=2, linestyle=\"--\")\n",
    "    ax.axvline(p_fit[4]+(p_fit[5]*cutoff),c=orange, linewidth=2, linestyle=\"--\")\n",
    "    ax.set_title(\"Line profile\",loc=\"left\",fontsize=14)\n",
    "    ax.set_xlabel(\"y (nm)\")\n",
    "    ax.set_ylabel(\"Counts\")  \n",
    "    ax.text(0.15,\n",
    "            0.7,\n",
    "            (\"Estimated Peaks:\\n\"\n",
    "            \"Peak 1 at {:.1f} nm\\n\"\n",
    "            \"Peak 2 at {:.1f} nm\\n\"\n",
    "            \"Fitted Peaks:\\n\"\n",
    "            \"Peak 1 at {:.1f} nm\\n\"\n",
    "            \"Peak 2 at {:.1f} nm\").format(peak1,peak2,p_fit[1],p_fit[4]),\n",
    "            horizontalalignment=\"center\",\n",
    "            verticalalignment=\"center\",\n",
    "            transform = ax.transAxes,\n",
    "            fontsize=12)\n",
    "\n",
    "    return ax\n",
    "\n",
    "\n",
    "def isolate_ring(data, c, w, axes=\"y\", cutoff=1.0):\n",
    "    \n",
    "    if axes == \"y\":\n",
    "        column = \"y_pick_rot\"\n",
    "    elif axes == \"x\":\n",
    "        column = \"x_pick_rot\"\n",
    "        \n",
    "    ring_data = data[np.where((data[column]>(c-(cutoff*w))) & (data[column] <(c+(cutoff*w))))]\n",
    "    return ring_data\n",
    "\n",
    "\n",
    "def rodrigues_rot(P, n0, n1):\n",
    "    # adapted from https://meshlogic.github.io/posts/jupyter/curve-fitting/fitting-a-circle-to-cluster-of-3d-points/\n",
    "    # If P is only 1d array (coords of single point), fix it to be matrix\n",
    "    if P.ndim == 1:\n",
    "        P = P[np.newaxis,:]\n",
    "    \n",
    "    # Get vector of rotation k and angle theta\n",
    "    n0 = n0/linalg.norm(n0)\n",
    "    n1 = n1/linalg.norm(n1)\n",
    "    k = np.cross(n0,n1)\n",
    "    k = k/linalg.norm(k)\n",
    "    theta =np.arccos(np.dot(n0,n1))\n",
    "    \n",
    "    # Compute rotated points\n",
    "    P_rot = np.zeros((len(P),3))\n",
    "    for i in range(len(P)):\n",
    "        P_rot[i] = P[i]*np.cos(theta) + np.cross(k,P[i])*np.sin(theta) + k*np.dot(k,P[i])*(1-np.cos(theta))\n",
    "\n",
    "    return P_rot\n",
    "\n",
    "\n",
    "def rotate_ring(XYZ): \n",
    "    # Fitting plane by SVD for the mean-centered data\n",
    "    # Eq. of plane is <p,n> + d = 0, where p is a point on plane and n is normal vector\n",
    "       \n",
    "    # Normal vector of fitting plane is given by 3rd column in V\n",
    "    # Note linalg.svd returns V^T, so we need to select 3rd row from V^T\n",
    "    ring_mean = XYZ.mean(axis=0)\n",
    "    ring_centered = XYZ - ring_mean\n",
    "    U,s,V = linalg.svd(ring_centered)\n",
    "    normal = V[2,:]\n",
    "    d = -np.dot(ring_mean, normal) \n",
    "        \n",
    "    n0 = normal # new z axes\n",
    "    n1 = [0,0,1] # old z axes\n",
    "\n",
    "    ring_rot = rodrigues_rot(ring_centered, n0, n1)\n",
    "    \n",
    "    return ring_rot, normal\n",
    "\n",
    "\n",
    "def plot_3d_ring(ring_rot, dist, ax=None):\n",
    "    \n",
    "    ax.scatter(ring_rot[0][:,0], ring_rot[0][:,1], ring_rot[0][:,2]-dist/2,c=blue,label=\"Spore\", alpha=0.5)\n",
    "    ax.scatter(ring_rot[1][:,0], ring_rot[1][:,1], ring_rot[1][:,2]+dist/2,c=orange,label=\"Mother\", alpha=0.5)\n",
    "    ax.set_xlabel(\"x (nm)\")\n",
    "    ax.set_ylabel(\"y (nm)\")\n",
    "    ax.set_zlabel(\"z (nm)\")\n",
    "    ax.legend(loc=\"best\",labelspacing=0.1)\n",
    "    #ax.view_init(elev=2., azim=10.)\n",
    "    ax.view_init(elev=30., azim=10.)\n",
    "    set_axes_equal_3d(ax)\n",
    "\n",
    "    return ax\n",
    "    \n",
    "    \n",
    "def set_axes_equal_3d(ax):\n",
    "    limits = np.array([ax.get_xlim3d(), ax.get_ylim3d(), ax.get_zlim3d()])\n",
    "    spans = abs(limits[:,0] - limits[:,1])\n",
    "    centers = np.mean(limits, axis=1)\n",
    "    radius = 0.5 * max(spans)\n",
    "    ax.set_xlim3d([centers[0]-radius, centers[0]+radius])\n",
    "    ax.set_ylim3d([centers[1]-radius, centers[1]+radius])\n",
    "    ax.set_zlim3d([centers[2]-radius, centers[2]+radius])\n",
    "    \n",
    "\n",
    "def angle_between(u, v, n=None):\n",
    "    if n is None:\n",
    "        return np.arctan2(np.linalg.norm(np.cross(u,v)), np.dot(u,v))*180/np.pi\n",
    "    else:\n",
    "        return np.arctan2(np.dot(n,np.cross(u,v)), np.dot(u,v))*180/np.pi\n",
    "\n",
    "\n",
    "def calc_R(x, y, xc, yc):\n",
    "    # adapted from https://gist.github.com/lorenzoriano/6799568\n",
    "    \"\"\" calculate the distance of each 2D points from the center (xc, yc) \"\"\"\n",
    "    return np.sqrt((x-xc)**2 + (y-yc)**2)\n",
    "\n",
    "\n",
    "def f(c, x, y):\n",
    "    \"\"\" calculate the algebraic distance between the data points and the mean circle centered at c=(xc, yc) \"\"\"\n",
    "    Ri = calc_R(x, y, *c)\n",
    "    return Ri - np.median(Ri)\n",
    "\n",
    "\n",
    "def leastsq_circle(x,y):\n",
    "    \n",
    "    x_m = np.median(x)\n",
    "    y_m = np.median(y)\n",
    "    center_estimate = x_m, y_m\n",
    "    center, ier = optimize.leastsq(f, center_estimate, args=(x,y))\n",
    "    xc, yc = center\n",
    "    Ri       = calc_R(x, y, *center)\n",
    "    R        = np.median(Ri)\n",
    "    residu   = np.sum((Ri - R)**2)\n",
    "    return xc, yc, R, residu\n",
    "\n",
    "\n",
    "def plot_data_circle(x, y, xc, yc, R, id, center=True, ax=None):\n",
    "    \n",
    "    if id == 0:\n",
    "        label = \"spore\"\n",
    "        color = blue\n",
    "    else:\n",
    "        label = \"mother\"\n",
    "        color = orange\n",
    "    \n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "    \n",
    "    if center:\n",
    "        x -= xc\n",
    "        y -= yc\n",
    "        xc = 0\n",
    "        yc = 0 \n",
    "        \n",
    "    #f, ax = plt.subplots(figsize=(5,5))  #figsize=(7, 5.4), dpi=72,\n",
    "    theta_fit = np.linspace(-np.pi,np.pi, 180)\n",
    "    x_fit = xc + R*np.cos(theta_fit)\n",
    "    y_fit = yc + R*np.sin(theta_fit)\n",
    "    \n",
    "    # plot fit\n",
    "    ax.plot(x_fit, y_fit, label=\"Fitted circle\", lw=2, c=red)\n",
    "    ax.plot([xc], [yc], mec=\"y\", mew=1,  c=red)\n",
    "    \n",
    "    # plot data\n",
    "    ax.scatter(x, y, alpha=0.4,  label=\"Projected locs\", marker=\".\", c=color)\n",
    "    \n",
    "    ax.set_xlabel(\"x rotated (nm)\")\n",
    "    ax.set_ylabel(\"y rotated (nm)\")\n",
    "    #ax.axis(\"equal\")\n",
    "    ax.set_xlim([-800,800])\n",
    "    ax.set_aspect(\"equal\",adjustable=\"datalim\")\n",
    "    \n",
    "    ax.legend(loc=\"best\", labelspacing=0.1)\n",
    "    ax.set_title(\"Least squares circle {}\\n\"\n",
    "                 \"Fit radius: {:.1f} nm\".format(label,R),loc=\"left\",fontsize=14)\n",
    "    \n",
    "    return \n",
    "\n",
    "\n",
    "def check_consistency(ring_par):\n",
    "    \n",
    "    ev = 0\n",
    "    \n",
    "    amp_1, c_1, w_1 = ring_par[0]\n",
    "    amp_2, c_2, w_2 = ring_par[1]\n",
    "    \n",
    "    if amp_1/amp_2 >= 0.3 and amp_2/amp_1 >= 0.3:\n",
    "        ev +=1 \n",
    "        \n",
    "    if abs(c_1-c_2) > 30 and abs(c_1-c_2) < 300:\n",
    "        ev +=1\n",
    "        \n",
    "    if w_1/w_2 >= 0.3 and w_2/w_1 >=0.3:\n",
    "        ev +=1\n",
    "        \n",
    "    return (ev==3)\n",
    "\n",
    "def plot_cum_exp(pooled_locs, fit_result_len, fit_result_dark, id, ax=None):\n",
    "    \n",
    "    if id == 0:\n",
    "        color = blue\n",
    "        label = 'spore'\n",
    "    if id == 1:\n",
    "        color = orange\n",
    "        label = 'mother'\n",
    "\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "    \n",
    "    data = pooled_locs.dark\n",
    "    data.sort()\n",
    "    y = np.arange(1, len(data) + 1)\n",
    "       \n",
    "    a = fit_result_dark.best_values[\"a\"]\n",
    "    t = fit_result_dark.best_values[\"t\"]\n",
    "    c = fit_result_dark.best_values[\"c\"]\n",
    "\n",
    "    ax.set_title(\n",
    "        \"Dark time (cumulative) {}\\n\"\n",
    "        r\"$Fit: {:.2f}\\cdot(1-exp(x/{:.2f}))+{:.2f}$\".format(label, a, t, c),loc=\"left\",fontsize=14)\n",
    "    data = pooled_locs.dark\n",
    "    data.sort()\n",
    "    y = np.arange(1, len(data) + 1)\n",
    "\n",
    "    ax.semilogx(data, y, c=color, label=\"Data\")\n",
    "    ax.semilogx(data, fit_result_dark.best_fit, c=red, label=\"Fit\")\n",
    "    ax.legend(loc=\"best\")\n",
    "    ax.set_xlabel(\"Duration (frames)\")\n",
    "    ax.set_ylabel(\"Frequency\")\n",
    "    \n",
    "    return ax\n",
    "\n",
    "\n",
    "def save_ring_locs(locs, info, path, file_id, pick, id, link=False):\n",
    "    \n",
    "    if link:\n",
    "        ending = \"_link.hdf5\"\n",
    "    else:\n",
    "        ending = \".hdf5\"\n",
    "       \n",
    "    locs.x /= 130\n",
    "    locs.y /= 130\n",
    "    \n",
    "    locs_name = \"file_{}_pick_{}_ring_{}{}\".format(file_id, pick, id, ending)\n",
    "    locs_path = os.path.join(path,\"ring_locs\")\n",
    "    locs_path_name = os.path.join(locs_path, locs_name)\n",
    "    \n",
    "    if not os.path.isdir(locs_path):\n",
    "        os.makedirs(locs_path)\n",
    "    \n",
    "    io.save_locs(locs_path_name, locs, info)\n",
    "\n",
    "    \n",
    "def export_pick_img(locs, path, file_id, pick, id, link=False):\n",
    "    \n",
    "    if link:\n",
    "        ending = \"_link.png\"\n",
    "    else:\n",
    "        ending = \".png\"\n",
    "    \n",
    "    pixelsize = 130\n",
    "\n",
    "    export_locs = locs.copy()\n",
    "    \n",
    "    export_locs.x /= pixelsize\n",
    "    export_locs.y /= pixelsize\n",
    "    \n",
    "    x_min = np.min(export_locs.x)    \n",
    "    x_max = np.max(export_locs.x)\n",
    "    y_min = np.min(export_locs.y)\n",
    "    y_max = np.max(export_locs.y)\n",
    "\n",
    "    viewport =  (y_min, x_min), (y_max, x_max)\n",
    "    oversampling = 50\n",
    "    len_x, image = render(export_locs, viewport = viewport, oversampling=oversampling, blur_method=\"smooth\")\n",
    "    \n",
    "    img_name = \"file_{}_pick_{}_ring_{}{}\".format(file_id, pick, id, ending)\n",
    "    img_path = os.path.join(path,\"ring_images\")\n",
    "    img_path_name = os.path.join(img_path,img_name)\n",
    "    \n",
    "    if not os.path.isdir(img_path):\n",
    "        os.makedirs(img_path)\n",
    "    \n",
    "    plt.imsave(img_path_name, image, cmap=\"hot\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7 HDF5 files found.\n",
      "7 HDF5 files found.\n"
     ]
    }
   ],
   "source": [
    "path = gui_fname()\n",
    "div_path = os.path.join(path, \"div\")\n",
    "fts_path = os.path.join(path, \"fts\")\n",
    "div_filenames = load_files(div_path)\n",
    "fts_filenames = load_files(fts_path)\n",
    "\n",
    "plotting = True\n",
    "max_dist = 130 #nanometer\n",
    "max_dark_time = 15 #frames\n",
    "binning = 50 # binning for peak histogram\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Main loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5692659921b64f72984a1affa941e268",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='Processing files', max=7.0, style=ProgressStyle(descripti…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5adf99acb75447289dd57363df001448",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='Processing picks', max=10.0, style=ProgressStyle(descript…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Miniconda3\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: divide by zero encountered in true_divide\n",
      "C:\\Miniconda3\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "101ff039c4db4453850f6c078fded7f7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='Processing picks', max=12.0, style=ProgressStyle(descript…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "ename": "IndexError",
     "evalue": "index 1 is out of bounds for axis 0 with size 1",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-3-462ba23b4da7>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     57\u001b[0m         \u001b[1;31m################\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     58\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 59\u001b[1;33m         \u001b[0mestimated_peaks\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mr_par\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mhist_data\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfind_peaks\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mpick_locs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbinning\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mbinning\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxes\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m\"y\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     60\u001b[0m         \u001b[0mfts_estimated_peak\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfts_r_par\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfts_hist_data\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfind_peak\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfts_pick_locs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbinning\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mbinning\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxes\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m\"y\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     61\u001b[0m         \u001b[1;31m# estimated peaks form MeanShift analysis to initialize guassian fits\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-1-9adfa049ae9b>\u001b[0m in \u001b[0;36mfind_peaks\u001b[1;34m(data, binning, axes)\u001b[0m\n\u001b[0;32m    117\u001b[0m     \u001b[0mpeaks\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msort\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mms\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcluster_centers_\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mNone\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;31m# assuming that there are two large peaks from the two rings\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    118\u001b[0m     \u001b[0mpeak1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfloat\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mpeaks\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[0m\n\u001b[1;32m--> 119\u001b[1;33m     \u001b[0mpeak2\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfloat\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mpeaks\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[0m\n\u001b[0m\u001b[0;32m    120\u001b[0m     \u001b[0mestimated_peaks\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m{\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[0mpeak1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[0mpeak2\u001b[0m\u001b[1;33m}\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    121\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mIndexError\u001b[0m: index 1 is out of bounds for axis 0 with size 1"
     ]
    }
   ],
   "source": [
    "# initialize containers\n",
    "ring_data, fts_ring_data = [], []\n",
    "ring_locs,fts_ring_locs = {}, {}\n",
    "ring_locs_linked, fts_ring_locs_linked = {}, {}\n",
    "ring_kinetics, fts_ring_kinetics = {}, {}\n",
    "ring_kinetics_fit, fts_ring_kinetics_fit = {}, {}\n",
    "ring_radii, fts_ring_radii = {}, {}\n",
    "ring_rot, fts_ring_rot = {}, {}\n",
    "ring_angles,fts_ring_angles = {}, {}\n",
    "center, fts_center = {}, {}\n",
    "\n",
    "rings_excluded, fts_rings_excluded = [], []\n",
    "ring_n_events, fts_ring_n_events = {}, {}\n",
    "circle_plots, fts_circle_plots = {}, {}\n",
    "file_id = 0\n",
    "\n",
    "analysis_folder = os.path.join(path, \"analysis\")\n",
    "analysis_folder_excluded = os.path.join(analysis_folder, \"excluded\")\n",
    "\n",
    "# image export settings\n",
    "img_format = \".png\"\n",
    "dpi = 100\n",
    "\n",
    "# cutoff is mutliplied to the sigma of the peak fit for the seperation of the two rings\n",
    "# smaller cutoff means smaller ring sections\n",
    "cutoff = 1.25\n",
    "\n",
    "#prepare analysis folder\n",
    "if not os.path.isdir(analysis_folder):\n",
    "    os.makedirs(analysis_folder)\n",
    "    os.makedirs(analysis_folder_excluded)\n",
    "\n",
    "#iterate over locs in directory:\n",
    "for filename in tqdm(div_filenames, desc=\"Processing files\"):\n",
    "    \n",
    "    file_id +=1\n",
    "    \n",
    "    #load DIV4A locs and convert distances from px to nm (Attention!)\n",
    "    locs, info, pixelsize = load_data(os.path.join(div_path, filename))\n",
    "    \n",
    "    #look for corresponding FtsZ file\n",
    "    fts_filename = filename.replace('DivIVA','ftsz')\n",
    "    try:\n",
    "        fts_locs, fts_info, fts_pixelsize = load_data(os.path.join(fts_path, fts_filename))\n",
    "    except:\n",
    "        print('ERROR: Could not find FtsZ file for {}\\n Check file names again.'.format(filename))\n",
    "    \n",
    "    # iterate over picks in a file\n",
    "    for pick in tqdm(np.unique(locs.group), desc=\"Processing picks\"):\n",
    "     \n",
    "        # select locs from pick\n",
    "        pick_locs = locs[locs.group == pick]\n",
    "        fts_pick_locs = fts_locs[fts_locs.group == pick]\n",
    "        ################\n",
    "        # (1) estimate the postition of the two rings using a histogram along the pick direction\n",
    "        #  a ring should yield a gaussian peak (2D projection in XY)\n",
    "        ################\n",
    "        \n",
    "        estimated_peaks, r_par, hist_data = find_peaks(pick_locs, binning=binning, axes=\"y\")\n",
    "        fts_estimated_peak, fts_r_par, fts_hist_data = find_peak(fts_pick_locs, binning=binning, axes=\"y\")\n",
    "        # estimated peaks form MeanShift analysis to initialize guassian fits\n",
    "        # r_par yields an array containing: [amplitued_1, center_1, width_1, amplitude_2, center_2, width_2]\n",
    "        \n",
    "        # seperate ring parameter to the corresponding ring\n",
    "        ring_parameter = []\n",
    "        for i in range(0, 4, 3):\n",
    "            ring_parameter.append(r_par[i: i+3])\n",
    "        \n",
    "        fts_ring_parameter = [] # to work as the double band script\n",
    "        fts_ring_parameter.append(fts_r_par)\n",
    "            \n",
    "        if check_consistency(ring_parameter):\n",
    "        # check three criteria:\n",
    "        # i   distance between the peaks (30 nm < distance < 300 nm)\n",
    "        # ii  width of the gaussian fits of the peaks (ratio needs to be larger then 0.2 in both cases)\n",
    "        # iii \"signal strength\" of both rings should be similar (amplitude of the gaussian peak should be larger then 0.3)\n",
    "            \n",
    "            # set up plot with gridspec\n",
    "            fig = plt.figure(figsize=(14, 16), constrained_layout=True)\n",
    "            gs = fig.add_gridspec(3, 2)\n",
    "            fig.suptitle((\"Pick {} - Ring analysis\\n\"\n",
    "                          \"File: {}\").format(pick,filename),\n",
    "                          fontsize=16,\n",
    "                          ha=\"center\")\n",
    "            \n",
    "            ax1 = fig.add_subplot(gs[0, 0])\n",
    "           \n",
    "            plot_peak_dist(pick_locs,\n",
    "                           hist_data,\n",
    "                           estimated_peaks,\n",
    "                           r_par,\n",
    "                           binning=binning,\n",
    "                           cutoff=cutoff,\n",
    "                           axes=\"y\",\n",
    "                           ax=ax1)\n",
    "            \n",
    "            peak_dist = abs(ring_parameter[0][1]-ring_parameter[1][1])\n",
    "            fts_peak_dist = abs(ring_parameter[0][1]-fts_ring_parameter[0][1])\n",
    "        \n",
    "            ################\n",
    "            # (2) isolate the locs from every ring\n",
    "            # fitted center +- peak-to-peak distance/2\n",
    "            ################\n",
    "            \n",
    "            #DIV4A\n",
    "            for id, ring in enumerate(ring_parameter):\n",
    "                \n",
    "                # isolate the locs from every ring\n",
    "                r_locs = isolate_ring(pick_locs, c=ring[1], w=ring[2], cutoff=cutoff, axes=\"y\")\n",
    "                ring_locs[id] = r_locs\n",
    "                \n",
    "                # link the isolated locs\n",
    "                r_locs_linked = link(r_locs, info, max_dist, max_dark_time)\n",
    "                ring_locs_linked[id] = r_locs_linked\n",
    "                \n",
    "                # extract number of locs\n",
    "                ring_n_events[id] = len(r_locs_linked)\n",
    "             \n",
    "            #FtsZ\n",
    "            for fts_id, fts_ring in enumerate(fts_ring_parameter):\n",
    "    \n",
    "                # isolate the locs from every ring\n",
    "                fts_r_locs = isolate_ring(fts_pick_locs, c=fts_ring[1], w=fts_ring[2], cutoff=cutoff, axes=\"y\")\n",
    "                fts_ring_locs[fts_id] = fts_r_locs\n",
    "                \n",
    "                # link the isolated locs\n",
    "                fts_r_locs_linked = link(fts_r_locs, fts_info, max_dist, max_dark_time)\n",
    "                fts_ring_locs_linked[fts_id] = fts_r_locs_linked\n",
    "                \n",
    "                # extract number of locs\n",
    "                fts_ring_n_events[fts_id] = len(fts_r_locs_linked)\n",
    "            \n",
    "              \n",
    "            ################\n",
    "            # (3b) calculate ring radius by rotation into plane of the ring and circle fitting in 2D projection\n",
    "            ################\n",
    "            \n",
    "            ax3 = fig.add_subplot(gs[1, 0])\n",
    "            ax5 = fig.add_subplot(gs[2, 0])\n",
    "            normal = {}\n",
    "            fts_normal = {}\n",
    "            \n",
    "            for id, rings in ring_locs.items():\n",
    "                \n",
    "                # prepare XYZ coordinates\n",
    "                XYZ = np.array([rings.x, rings.y, rings.z]).transpose()\n",
    "\n",
    "                # rotate coordinates into the plane of the ring. ring should be now in X-Y\n",
    "                ring_rot[id], normal[id] = rotate_ring(XYZ)\n",
    "                \n",
    "                # fit 2D ring into the XY ring data\n",
    "                x_center, y_center, Radius, residu = leastsq_circle(ring_rot[id][:,0], ring_rot[id][:,1])\n",
    "                \n",
    "                c = rodrigues_rot(np.array([x_center, y_center, 0]), [0,0,1], normal[id]) + XYZ.mean(axis=0)\n",
    "                center[id] = c.flatten()\n",
    "                \n",
    "                # extract radii\n",
    "                ring_radii[id] = Radius\n",
    "                \n",
    "                # plot\n",
    "                if id == 0:\n",
    "                    ax_circ = ax3\n",
    "                if id == 1:\n",
    "                    ax_circ = ax5\n",
    "                \n",
    "                plot_data_circle(ring_rot[id][:,0],\n",
    "                                 ring_rot[id][:,1],\n",
    "                                 x_center,\n",
    "                                 y_center,\n",
    "                                 Radius,\n",
    "                                 id,\n",
    "                                 center=True,\n",
    "                                 ax=ax_circ)\n",
    "                \n",
    "            for fts_id, fts_rings in fts_ring_locs.items():\n",
    "                \n",
    "                # prepare XYZ coordinates\n",
    "                fts_XYZ = np.array([fts_rings.x, fts_rings.y, fts_rings.z]).transpose()\n",
    "\n",
    "                # rotate coordinates into the plane of the ring. ring should be now in X-Y\n",
    "                fts_ring_rot[fts_id], fts_normal[fts_id] = rotate_ring(fts_XYZ)\n",
    "                \n",
    "                # fit 2D ring into the XY ring data\n",
    "                fts_x_center, fts_y_center, fts_Radius, fts_residu = leastsq_circle(fts_ring_rot[fts_id][:,0],\n",
    "                                                                                    fts_ring_rot[fts_id][:,1])\n",
    "                \n",
    "                fts_c = rodrigues_rot(np.array([fts_x_center, fts_y_center,0]), \n",
    "                                      [0,0,1], \n",
    "                                      fts_normal[fts_id]) + fts_XYZ.mean(axis=0)\n",
    "                \n",
    "                fts_center[fts_id] = fts_c.flatten()\n",
    "                \n",
    "                # extract radii\n",
    "                fts_ring_radii[fts_id] = fts_Radius\n",
    "                \n",
    "            \n",
    "            # calculate the center-to-center distance\n",
    "            center_dist = distance.euclidean(center[0], center[1])\n",
    "            fts_center_dist = distance.euclidean(center[0], fts_center[0])\n",
    "                \n",
    "            # calculate the angle between to two fitted ring planes\n",
    "            ring_angle = abs(angle_between(normal[0],normal[1]))\n",
    "            \n",
    "            # if normals point into opposite direction\n",
    "            if ring_angle > 90:\n",
    "                ring_angle -=180\n",
    "                ring_angle = abs(ring_angle)\n",
    "               \n",
    "            # calculate the angle between coverslip and spore or mother ring\n",
    "            coverslip_s_angle = angle_between(normal[0],[0,0,1])\n",
    "            coverslip_m_angle = angle_between(normal[1],[0,0,1])\n",
    "            \n",
    "            # plot 3D ring\n",
    "            ax2 = fig.add_subplot(gs[0, 1], projection=\"3d\")\n",
    "            plot_3d_ring(ring_rot,center_dist, ax=ax2)\n",
    "            \n",
    "            ################\n",
    "            # (4) qPAINT analysis\n",
    "            ################\n",
    "            \n",
    "            ax4 = fig.add_subplot(gs[1, 1])\n",
    "            ax6 = fig.add_subplot(gs[2, 1])\n",
    "            \n",
    "            for id, rings_linked in ring_locs_linked.items():\n",
    "\n",
    "                    # compute kinetics\n",
    "                    kinetics = compute_dark_times(rings_linked)\n",
    "\n",
    "                    # compute mean bright and dark times\n",
    "                    ring_kinetics[id] = {\"bright\":np.nanmean(kinetics.len),\n",
    "                                         \"dark\":np.nanmean(kinetics.dark)}\n",
    "                    \n",
    "                    fit_result_len = fit_cum_exp(kinetics.len)\n",
    "                    fit_result_dark = fit_cum_exp(kinetics.dark)\n",
    "                \n",
    "                    if id == 0:\n",
    "                        ax_fit=ax4\n",
    "                    if id == 1:\n",
    "                        ax_fit=ax6\n",
    "                    \n",
    "                    plot_cum_exp(kinetics, fit_result_len, fit_result_dark, id, ax=ax_fit)\n",
    "                    \n",
    "                    ring_kinetics_fit[id] = {\"bright\":np.nanmean(fit_result_len.best_values[\"t\"]),\n",
    "                                             \"dark\":np.nanmean(fit_result_dark.best_values[\"t\"])}\n",
    "                    \n",
    "            for fts_id, fts_rings_linked in fts_ring_locs_linked.items():\n",
    "\n",
    "                    # compute kinetics\n",
    "                    fts_kinetics = compute_dark_times(fts_rings_linked)\n",
    "\n",
    "                    # compute mean bright and dark times\n",
    "                    fts_ring_kinetics[fts_id] = {\"bright\":np.nanmean(fts_kinetics.len),\n",
    "                                         \"dark\":np.nanmean(fts_kinetics.dark)}\n",
    "                    \n",
    "                    fts_fit_result_len = fit_cum_exp(fts_kinetics.len)\n",
    "                    fts_fit_result_dark = fit_cum_exp(fts_kinetics.dark)\n",
    "                    \n",
    "                    fts_ring_kinetics_fit[fts_id] = {\"bright\":np.nanmean(fts_fit_result_len.best_values[\"t\"]),\n",
    "                                             \"dark\":np.nanmean(fts_fit_result_dark.best_values[\"t\"])}\n",
    "                       \n",
    "            ################\n",
    "            # (5) append data into large array\n",
    "            ################\n",
    "            \n",
    "            ring_data.append([file_id, #running file index\n",
    "                              filename, #DIV4A filename,\n",
    "                              fts_filename, #FtsZ filename\n",
    "                              pick, #pick number\n",
    "                              ring_radii[0], #fitted radius of spore (nm)\n",
    "                              ring_radii[1], #fitted radius of mother (nm)\n",
    "                              fts_ring_radii[0], #fitted radius of ftsZ (nm)\n",
    "                              peak_dist, #peak-to-peak distance (nm)\n",
    "                              fts_peak_dist, #spore_peak-to-ftsZ_peak (nm)\n",
    "                              center_dist, #center-to-center distance (nm)\n",
    "                              fts_center_dist, #spore_center-to-ftsZ_center (nm)\n",
    "                              ring_angle, #angle between the two ring planes (°)\n",
    "                              coverslip_s_angle, #anle between spore ring and coverslip (°)\n",
    "                              coverslip_m_angle, #anle between mother ring and coverslip (°)\n",
    "                              ring_n_events[0], #number of locs in spore ring\n",
    "                              ring_n_events[1], #number of  locs in mother ring\n",
    "                              fts_ring_n_events[0], #number of locs in ftsZ ring\n",
    "                              ring_kinetics[0][\"bright\"], #mean bright time of spore ring (frames)\n",
    "                              ring_kinetics[0][\"dark\"], #mean dark time of spore ring (frames)\n",
    "                              ring_kinetics[1][\"bright\"], #mean bright time of mother ring (frames)\n",
    "                              ring_kinetics[1][\"dark\"], #mean dark time of mother ring (frames)\n",
    "                              fts_ring_kinetics[0][\"bright\"], #mean bright time of ftsZ ring (frames)\n",
    "                              fts_ring_kinetics[0][\"dark\"], #mean dark time of ftsZ ring (frames)\n",
    "                              ring_kinetics_fit[0][\"bright\"], #cumulative mean bright times of spore (frames)\n",
    "                              ring_kinetics_fit[0][\"dark\"], #cumulative mean dark times of spore (frames)\n",
    "                              ring_kinetics_fit[1][\"bright\"], #cumulative mean bright times of mother (frames)\n",
    "                              ring_kinetics_fit[1][\"dark\"], #cumulative mean dark times of mother (frames)\n",
    "                              fts_ring_kinetics_fit[0][\"bright\"], #cumulative mean bright times of ftsZ (frames)\n",
    "                              fts_ring_kinetics_fit[0][\"dark\"], #cumulative mean dark times of ftsZ (frames)\n",
    "                             ])\n",
    "            \n",
    "            ################\n",
    "            # (6) save plots\n",
    "            ################\n",
    "            \n",
    "            img_fname = \"file_{}_pick_{}\".format(file_id, pick)\n",
    "            img_name = os.path.join(analysis_folder, img_fname)\n",
    "            plt.savefig(img_name+img_format, dpi=dpi, format=\"png\")\n",
    "            plt.close(fig)\n",
    "            \n",
    "            \n",
    "            ################\n",
    "            # (7) export ring locs\n",
    "            ################\n",
    "            \n",
    "            for id, rings_linked in ring_locs_linked.items():\n",
    "                save_ring_locs(rings_linked, info, analysis_folder, file_id, pick, id, link=True)\n",
    "            \n",
    "            for id, rings in ring_locs.items():\n",
    "                save_ring_locs(rings, info, analysis_folder, file_id, pick, id, link=False)\n",
    "            \n",
    "        else:\n",
    "            fig, ax = plt.subplots(1, 1, figsize=(7, 5), constrained_layout=True)\n",
    "            fig.suptitle((\"Pick {} - Ring analysis\\n\"\n",
    "                          \"File: {}\").format(pick,filename), fontsize=14, ha=\"center\")\n",
    "            \n",
    "            plot_peak_dist(pick_locs, hist_data, estimated_peaks, r_par, binning=binning, axes=\"y\", ax=ax)\n",
    "            img_fname = \"file_{}_pick_{}\".format(file_id, pick)\n",
    "            img_name = os.path.join(analysis_folder_excluded, img_fname)\n",
    "            plt.savefig(img_name+img_format, format=\"png\")\n",
    "            plt.close(fig)\n",
    "            rings_excluded.append(pick)\n",
    "        \n",
    "print(\"Calculation finished.\")\n",
    "print(\"Total: {} rings analysed, {} excluded.\".format(len(ring_data), len(rings_excluded)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save analysis data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_ring_data = pd.DataFrame(ring_data, columns=[\"file_index\",\n",
    "                                                \"div filename\",\n",
    "                                                \"fts filename\",\n",
    "                                                \"group\",\n",
    "                                                \"radius_spore\",\n",
    "                                                \"radius_mother\",\n",
    "                                                \"radius_ftsZ\",\n",
    "                                                \"p2p distance\",\n",
    "                                                \"p2p spore-ftsZ distance\",\n",
    "                                                \"c2c distance\",\n",
    "                                                \"c2c spore-ftsZ distance\",\n",
    "                                                \"angle_between_rings\",\n",
    "                                                \"angle_between_s_ring_and_coverslip\",\n",
    "                                                \"angle_between_m_ring_and_coverslip\",\n",
    "                                                \"n_events_spore\",\n",
    "                                                \"n_events_mother\",\n",
    "                                                \"n_events_ftsZ\",\n",
    "                                                \"mean_bright_spore\",\n",
    "                                                \"mean_dark_spore\",\n",
    "                                                \"mean_bright_mother\",\n",
    "                                                \"mean_dark_mother\",\n",
    "                                                \"mean_bright_ftsZ\",\n",
    "                                                \"mean_dark_ftsZ\",\n",
    "                                                \"fit_bright_spore\",\n",
    "                                                \"fit_dark_spore\",\n",
    "                                                \"fit_bright_mother\",\n",
    "                                                \"fit_dark_mother\",\n",
    "                                                \"fit_bright_ftsZ\",\n",
    "                                                \"fit_dark_ftsZ\"\n",
    "                                               ])\n",
    "df_ring_data.to_csv(\"exchange_data.csv\")\n",
    "print(\"Data saved to CSV file in locs folder.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "spor-paint",
   "language": "python",
   "name": "spor-paint"
  },
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}