dual-band.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 = float(peaks[0])\n",
    "    peak2 = 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",
    "\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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = gui_fname()\n",
    "filenames = load_files(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": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4db0dd2e46dd4c63a1d2224243e3a2c4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing files:   0%|          | 0/61 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ed207219345f4af4a98bf0743c507694",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/14 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: divide by zero encountered in true_divide\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: invalid value encountered in true_divide\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: divide by zero encountered in true_divide\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c57f05c8d1514dc3b390e43d370b3f4f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/14 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\scipy\\optimize\\minpack.py:834: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  category=OptimizeWarning)\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: divide by zero encountered in true_divide\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ba28c636b9e7483f8d3ae8aa1abe7ad3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/10 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8674841c70f2409c90a222f219e29e60",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/17 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\scipy\\optimize\\minpack.py:834: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  category=OptimizeWarning)\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4b35f425bf18414b8408c9497294888c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/17 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: divide by zero encountered in true_divide\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: invalid value encountered in true_divide\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: divide by zero encountered in true_divide\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "61337941b3714bf28dfa952b08534705",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/13 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\scipy\\optimize\\minpack.py:834: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  category=OptimizeWarning)\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0f5151385b83426ebef34fa44c812124",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/15 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d6e099b98beb45f98fdfca705c33dc16",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/19 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: divide by zero encountered in true_divide\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "64b6ba704b6e4a11afdf1bb879b07c61",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/15 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: divide by zero encountered in true_divide\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3bb7604bb78a4f88843934c5a20e6eb6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/1 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9ee7a852c64f4d5db77a8b093dfdfc94",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/3 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "784444f9e5f8411f86acd1f3a1ce2d8d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/1 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "81ffabd8eff94aa49ed48b37f48cda82",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/4 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b86b8c5b193b40e5afa0ebb9af000771",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/1 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "46e3810545874fc79dba0e5a4b439cc3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/3 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: divide by zero encountered in true_divide\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: invalid value encountered in true_divide\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: divide by zero encountered in true_divide\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2f91e769de294aa1b26c147c4504bcbf",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/1 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: divide by zero encountered in true_divide\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3bb3eb3dec104e2bb87718a3ef5ffcb8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/2 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "51331cfe9d5f4892b4c11f90d500a3de",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/1 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: divide by zero encountered in true_divide\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1c32b9c388904907b3afb8004b6b46c5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/1 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "198b9f661f8944368dfa166f8c3320d5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/5 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: divide by zero encountered in true_divide\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "faf9087aede74c4cb667a9a76c77cba2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/1 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b2073ec8e5a547cf8129806aa6b73c12",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/6 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4d23bce05fb14727b40e358cb626f82c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/3 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bd1b532f33a94a8f84c76900ec51821d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/7 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\scipy\\optimize\\minpack.py:834: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  category=OptimizeWarning)\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: divide by zero encountered in true_divide\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: invalid value encountered in true_divide\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: divide by zero encountered in true_divide\n",
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\ipykernel_launcher.py:99: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "08b254c611d042f5938e9f275f3647df",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/5 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\data\\reinhardt\\envs\\spor-paint\\lib\\site-packages\\scipy\\optimize\\minpack.py:834: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  category=OptimizeWarning)\n"
     ]
    },
    {
     "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~\\AppData\\Local\\Temp\\2\\ipykernel_29932\\1719932740.py\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     50\u001b[0m         \u001b[1;31m################\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     51\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 52\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[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     53\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[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     54\u001b[0m         \u001b[1;31m# r_par yields an array containing: [amplitued_1, center_1, width_1, amplitude_2, center_2, width_2]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\AppData\\Local\\Temp\\2\\ipykernel_29932\\798507837.py\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[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    118\u001b[0m     \u001b[0mpeak1\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[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 119\u001b[1;33m     \u001b[0mpeak2\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[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[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 = []\n",
    "ring_locs = {}\n",
    "ring_locs_linked = {}\n",
    "ring_kinetics = {}\n",
    "ring_kinetics_fit = {}\n",
    "ring_radii = {}\n",
    "ring_rot = {}\n",
    "ring_angles = {}\n",
    "center = {}\n",
    "\n",
    "rings_excluded = []\n",
    "ring_n_events = {}\n",
    "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(filenames, desc=\"Processing files\"):\n",
    "    \n",
    "    file_id +=1\n",
    "    \n",
    "    #load locs and convert distances from px to nm (Attention!)\n",
    "    locs, info, pixelsize = load_data(os.path.join(path,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",
    "\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",
    "        # 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",
    "            \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",
    "        \n",
    "            ################\n",
    "            # (2) isolate the locs from every ring\n",
    "            # fitted center +- peak-to-peak distance/2\n",
    "            ################\n",
    "            \n",
    "            \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",
    "                        \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",
    "            for id, rings in ring_locs.items():\n",
    "                \n",
    "                # prepare XYZ coordinates \n",
    "                # here x,y,z coordinates are the default locs coordinates, not the rotated picked locs 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",
    "                # rotate centers back into default coodinate system\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",
    "            # calculate the center-to-center distance\n",
    "            center_dist = distance.euclidean(center[0], center[1])\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",
    "            ################\n",
    "            # (5) append data into large array\n",
    "            ################\n",
    "            \n",
    "            ring_data.append([file_id, #running file index\n",
    "                              filename, #filename\n",
    "                              pick, #pick number\n",
    "                              ring_radii[0], #fitted radius of spore (nm)\n",
    "                              ring_radii[1], #fitted radius of mother (nm)\n",
    "                              peak_dist, #peak-to-peak distance (nm)\n",
    "                              center_dist, #center-to-center distance (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",
    "                              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",
    "                              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",
    "                             ])\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": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data saved to CSV file in locs folder.\n"
     ]
    }
   ],
   "source": [
    "df_ring_data = pd.DataFrame(ring_data, columns=[\"file_index\",\n",
    "                                                \"filename\",\n",
    "                                                \"group\",\n",
    "                                                \"radius_spore\",\n",
    "                                                \"radius_mother\",\n",
    "                                                \"p2p distance\",\n",
    "                                                \"c2c 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",
    "                                                \"mean_bright_spore\",\n",
    "                                                \"mean_dark_spore\",\n",
    "                                                \"mean_bright_mother\",\n",
    "                                                \"mean_dark_mother\",\n",
    "                                                \"fit_bright_spore\",\n",
    "                                                \"fit_dark_spore\",\n",
    "                                                \"fit_bright_mother\",\n",
    "                                                \"fit_dark_mother\"\n",
    "                                               ])\n",
    "df_ring_data.to_csv(\"ring_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": "Python 3 (ipykernel)",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}