exchange-paint-FtsZ-ZapA.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",
    "import re\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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 load_ring_data_df(dirname, filename):\n",
    "    file = ospath.join(dirname, filename)\n",
    "\n",
    "    try: \n",
    "        df = pd.read_pickle(file)\n",
    "    except FileNotFoundError:\n",
    "        print(\"No results of previously analyzed datasets were detected.\")\n",
    "        return None\n",
    "    else: \n",
    "        print(\"Results of previously analyzed datasets were detected.\")\n",
    "        return df\n",
    "        \n",
    "\n",
    "def identify_new_files(files, df_ring_data):\n",
    "    \"\"\"\n",
    "    Identify which files have already been analyzed previously. \n",
    "    Return list of new files for processing\n",
    "    \"\"\"\n",
    "    new_files = []\n",
    "    for file in files:\n",
    "        if file not in df_ring_data['filename'].values:\n",
    "            new_files.append(file)\n",
    "    \n",
    "    n_old = len(files)-len(new_files)\n",
    "    n_new = len(new_files)\n",
    "    if n_old == 1:\n",
    "        print(\" {} HDF5 file was previously analyzed.\".format(n_old))\n",
    "    else: \n",
    "        print(\" {} HDF5 files were previously analyzed.\".format(n_old))\n",
    "        \n",
    "    if n_new == 1:\n",
    "        print(\" {} HDF5 file is new and will be analyzed.\".format(n_new))\n",
    "    else: \n",
    "        print(\" {} HDF5 files are new and will be analyzed.\".format(n_new))\n",
    "        \n",
    "    return new_files\n",
    "   \n",
    "\n",
    "def identify_fov_cell_type_channel(df_ring_data, filenames, fov_id_start, protein1, protein2):\n",
    "    \"\"\"\n",
    "    For each file identify \n",
    "    - cell type: sporulating or vegetative cells\n",
    "    - fov index: fov from which the picks were generated\n",
    "    Results will be saved in a dataframe with columns:\n",
    "    'fov_id', 'filename', 'cell_type', 'channel'\n",
    "    \n",
    "    If some files have been analyzed before the current script execution \n",
    "    the results were saved in the ring_data file and loaded to df_ring_data.\n",
    "    (columns: \"fov_id\", \"cell_type\", \"filename\", \"group\", ... where \n",
    "    group are the pick ids.)\n",
    "    If no prior analysis results exist, then df_ring_data = None. \n",
    "    Therefore we can check if a new file belongs to a previously analyzed fov.\n",
    "    \n",
    "    \"\"\"\n",
    "    dictionary = {}\n",
    "    filenames_no_cell_type = []\n",
    "    filenames_no_channel = []\n",
    "\n",
    "    \n",
    "    fov_id_counter = fov_id_start + 1\n",
    "    \n",
    "    for i, filename in enumerate(filenames):\n",
    "        # cell type: spor or veg:\n",
    "        cell_type = np.nan\n",
    "                \n",
    "        spor_found = re.search('spor', filename, re.IGNORECASE)\n",
    "        veg_found = re.search('veg', filename, re.IGNORECASE)\n",
    "        \n",
    "        if spor_found and not veg_found:\n",
    "            cell_type = 'spor'\n",
    "        elif not spor_found and veg_found:\n",
    "            cell_type = 'veg'\n",
    "        elif spor_found and veg_found:\n",
    "            # consider the string occuring first as cell type determining string\n",
    "            spor_found_location = spor_found.start()\n",
    "            veg_found_location = veg_found.start()\n",
    "            \n",
    "            if spor_found_location < veg_found_location:\n",
    "                cell_type = 'spor'\n",
    "            else:\n",
    "                cell_type = 'veg'\n",
    "        else:\n",
    "            cell_type = np.nan\n",
    "            filenames_no_cell_type.append(filename)\n",
    "            \n",
    "        \n",
    "        # channel: protein1 or protein2:\n",
    "        channel = np.nan\n",
    "                \n",
    "        protein1_found = re.search(protein1, filename, re.IGNORECASE)\n",
    "        protein2_found = re.search(protein2, filename, re.IGNORECASE)\n",
    "        \n",
    "        if protein1_found and not protein2_found:\n",
    "            channel = protein1\n",
    "        elif not protein1_found and protein2_found:\n",
    "            channel = protein2\n",
    "        elif protein1_found and protein2_found:\n",
    "            # consider the string occuring first as channel determining string\n",
    "            protein1_found_location = protein1_found.start()\n",
    "            protein2_found_location = protein2_found.start()\n",
    "            \n",
    "            if protein1_found_location < protein2_found_location:\n",
    "                channel = protein1\n",
    "            else:\n",
    "                channel = protein2\n",
    "        else:\n",
    "            channel = np.nan\n",
    "            filenames_no_channel.append(filename)\n",
    "\n",
    "        \n",
    "        # search if a file from the same fov was already registered:\n",
    "        if not pd.isnull(cell_type):\n",
    "            \"\"\"\n",
    "            The string starting after the second occurance of '_' and ending with 'kcb102' \n",
    "            is the same for all files from the same FOV, no meatter their cell type or \n",
    "            the protein channel.\n",
    "            However, this might be the same string for other FOVs measured on the same day. \n",
    "            Thus, also search for the string 'fov' in the filename, which indicates the FOV id of that specific day.\n",
    "            \"\"\"\n",
    "            \n",
    "            # get substrings after the second '_'\n",
    "            kcb102_idx = filename.find('kcb102')\n",
    "            len_kcb = len('kcb102')\n",
    "            if kcb102_idx == -1:\n",
    "                kcb102_idx = filename.find('kcb1113')\n",
    "                len_kcb = len('kcb1113')\n",
    "                if kcb102_idx == -1:\n",
    "                    kcb102_idx = filename.find('kcb328')\n",
    "                    len_kcb = len('kcb328')\n",
    "            fov_identifying_substring = filename[second_underscore_idx:kcb102_idx+len_kcb]\n",
    "            # Is there the word 'fov' in the filename?\n",
    "            # If yes, then extract the number after 'fov' to search other files with this name\n",
    "            if re.search('fov', filename, re.IGNORECASE):\n",
    "                fov_filename_id = re.findall(r'%s(\\d+)' % 'fov', filename.lower())\n",
    "                if len(fov_filename_id) == 1:\n",
    "                    fov_substring = 'fov' + str(fov_filename_id[0])\n",
    "                else:\n",
    "                    fov_substring = ''\n",
    "            #print(filename)\n",
    "            #print(fov_identifying_substring)\n",
    "            #print(fov_substring)\n",
    "            \n",
    "            \n",
    "            # check if an already registered file exists that contains the substrings\n",
    "            filenames_found = []\n",
    "            old_or_new_file = [] # True if file from previous run of the script, False if new file.\n",
    "            for filename2 in dictionary.keys(): # search in files that were already registered.\n",
    "                #print('  ', filename2)\n",
    "                contained = fov_identifying_substring.lower() in filename2.lower() and fov_substring.lower() in filename2.lower()\n",
    "                #print('  ', fov_identifying_substring.lower() in filename2.lower())\n",
    "                #print('  ', fov_substring.lower() in filename2.lower())\n",
    "                if contained and filename != filename2:\n",
    "                    filenames_found.append(filename2)\n",
    "                    old_or_new_file.append(False)\n",
    "            if isinstance(df_ring_data, pd.DataFrame):\n",
    "                previously_analyzed_files = np.unique(df_ring_data['filename'])\n",
    "                for filename2 in previously_analyzed_files:\n",
    "                    contained = fov_identifying_substring.lower() in filename2.lower() and fov_substring.lower() in filename2.lower()\n",
    "                    if contained and filename != filename2:\n",
    "                        filenames_found.append(filename2)\n",
    "                        old_or_new_file.append(True)\n",
    "            #print()\n",
    "            \n",
    "            # if one other file was found: assign the existing fov index to the newly registred file\n",
    "            # if no other file was found: assign a new fov index to the newly registered file\n",
    "            if len(filenames_found) > 3:\n",
    "                raise Exception('''Files from the same FOV than ''' + filename + ''' where searched. \n",
    "                However more than three other file was detected: \n",
    "                ''' + '\\n'.join(filenames_found))\n",
    "            elif len(filenames_found) >= 1:\n",
    "                \n",
    "                filename_found = filenames_found[0]\n",
    "                if old_or_new_file[0]: # filename_found is from previous run of the script\n",
    "                    fov_id = df_ring_data.loc[df_ring_data['filename'] == filename_found, 'fov_id'].iloc[0]\n",
    "                if not old_or_new_file[0]: # filename_found is also a new file.\n",
    "                    fov_id = dictionary[filename_found][0]\n",
    "            else: # No file from the same fov previously registered\n",
    "                fov_id = fov_id_counter\n",
    "                fov_id_counter += 1\n",
    "        \n",
    "        else:\n",
    "            fov_id = np.nan\n",
    "        \n",
    "        \n",
    "        \n",
    "        dictionary[filename] = [fov_id, filename, cell_type, channel]\n",
    "        \n",
    "    df_results = pd.DataFrame.from_dict(dictionary, orient = 'index', columns = ['fov_id', 'filename', 'cell_type', 'channel'])\n",
    "    df_results = df_results.reset_index()\n",
    "\n",
    "    print()\n",
    "    print('The cell type (spr or veg) of these files could not be determined')\n",
    "    print('and thus cannot be used for further analysis:')\n",
    "    for filename in filenames_no_cell_type:\n",
    "        print(' -', filename)\n",
    "        \n",
    "    return df_results\n",
    "                \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_peak(data, binning=100, bins = None, 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 = float(ms.cluster_centers_[0])\n",
    "    estimated_peak = {0:peak1}\n",
    "\n",
    "    # histogram\n",
    "    if bins is not None:\n",
    "        n, bins = np.histogram(data[column], bins=bins)\n",
    "    else:\n",
    "        n, bins = np.histogram(data[column], bins=binning)\n",
    "    #bins = bins[1:]\n",
    "    centers = (bins[:-1] + bins[1:]) / 2\n",
    "    hist_data = [n, centers]\n",
    "\n",
    "    # fit peaks\n",
    "    p0 = [peak1/2, peak1, 40]\n",
    "    try:\n",
    "        p_fit, p_cov = optimize.curve_fit(gaus, centers, 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",
    "def plot_peak_dist(data, \n",
    "                   hist_data, \n",
    "                   peaks, \n",
    "                   p_fit, \n",
    "                   binning=100,\n",
    "                   cutoff=1, \n",
    "                   axes=\"y\", \n",
    "                   title = None, \n",
    "                   color_hist = gray, \n",
    "                   color_fit = red, \n",
    "                   annotate = True, \n",
    "                   ax=None):\n",
    "    \n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "        \n",
    "    \n",
    "    if axes == \"xyz\":\n",
    "        column == 1\n",
    "    \n",
    "    column = axes\n",
    "        \n",
    "    peak1 = peaks[0]\n",
    "\n",
    "    n = hist_data[0]\n",
    "    bins = hist_data[1]\n",
    "    \n",
    "    #fig, ax = plt.subplots(figsize=(9, 5))\n",
    "    binwidth = bins[1]-bins[0]\n",
    "    ax.bar(bins, n, width=binwidth, color=color_hist, alpha = 0.5)\n",
    "    xlin = np.linspace(0, data[column].max(), 1000)\n",
    "    ax.plot(xlin, gaus(xlin,*p_fit[0:3]), c=color_fit, 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",
    "    if title is None:\n",
    "        ax.set_title(\"Line profile\",loc=\"left\",fontsize=14)\n",
    "    else:\n",
    "        ax.set_title(title,loc=\"left\",fontsize=14)\n",
    "        \n",
    "    ax.set_xlabel(\"y (nm)\")\n",
    "    ax.set_ylabel(\"Counts\")  \n",
    "    if annotate:\n",
    "        ax.text(0.15,\n",
    "                0.7,\n",
    "                (\"Estimated Peaks:\\n\"\n",
    "                \"Peak at {:.1f} nm\\n\"\n",
    "                \"Fitted Peaks:\\n\"\n",
    "                \"Peak at {:.1f} nm\\n\"\n",
    "                \"sigma of {:.1f} nm\\n\").format(peak1, p_fit[1], p_fit[2]),\n",
    "                horizontalalignment=\"center\",\n",
    "                verticalalignment=\"center\",\n",
    "                transform = ax.transAxes,\n",
    "                fontsize=12)\n",
    "\n",
    "    return ax\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_rot1, ring_rot2, color1, color2, label1, label2, ax=None):\n",
    "    \n",
    "    ax.scatter(ring_rot1.x_pick_rot, ring_rot1.y_pick_rot, ring_rot1.z,c=color1, alpha=0.5, label = label1)\n",
    "    ax.scatter(ring_rot2.x_pick_rot, ring_rot2.y_pick_rot, ring_rot2.z,c=color2, alpha=0.5, label = label2)\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=30.)\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 (True)\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": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2 HDF5 files found.\n"
     ]
    }
   ],
   "source": [
    "#path = gui_fname()\n",
    "#path = r'X:\\users\\kcramer\\sporPAINT\\ZapA_FtsZ\\picked'\n",
    "#path = r'X:\\users\\kcramer\\sporPAINT\\SepF-FtsZ\\test_picked'\n",
    "path = r'X:\\users\\kcramer\\sporPAINT\\ZapA_FtsZ\\picked\\new_picks_6'\n",
    "\n",
    "protein1 = 'FtsZ'\n",
    "protein2 = 'ZapA'\n",
    "\n",
    "filenames_all = 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": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No results of previously analyzed datasets were detected.\n",
      "\n",
      "The cell type (spr or veg) of these files could not be determined\n",
      "and thus cannot be used for further analysis:\n"
     ]
    }
   ],
   "source": [
    "# Check if some of the found hdf5 files were already analyzed?\n",
    "# If yes, open ring_data dataframe with previous results.\n",
    "df_ring_data = load_ring_data_df(path, \"exchange_data.pkl\")\n",
    "\n",
    "\n",
    "# Identify which files have not yet been analyzed.\n",
    "if df_ring_data is not None:\n",
    "    filenames = identify_new_files(filenames_all, df_ring_data)\n",
    "    fov_id_start = df_ring_data['fov_id'].max()\n",
    "else:\n",
    "    filenames = filenames_all\n",
    "    fov_id_start = 0\n",
    "\n",
    "# Create a dataframe that saves which file was taken from which FOV and which cell types are contained (spor or veg)\n",
    "# and which channel the data is from (protein1 or protein2)\n",
    "# columns = ['fov_id', 'filename', 'cell_type', 'channel']\n",
    "df_fov_file_assign = identify_fov_cell_type_channel(df_ring_data, filenames, fov_id_start, protein1, protein2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FOV ID: 1\n",
      "  spor : \n",
      "    FtsZ : --\n",
      "    ZapA : --\n",
      "  veg : \n",
      "    FtsZ : veg_FtsZ_2plex_rings_230621_apollo_kcb328_fov1_2plex_r4-170pM_DP_FtsZ_1_drift_aligned_picked.hdf5\n",
      "    ZapA : veg_ZapA_2plex_rings_230621_apollo_kcb328_fov1_2plex_r3-1nM_DP_ZapA_1_drift_aligned_picked.hdf5\n"
     ]
    }
   ],
   "source": [
    "for fov_id in range(int(fov_id_start)+1, int(df_fov_file_assign['fov_id'].max())+1):\n",
    "    \n",
    "    print('FOV ID:', fov_id)\n",
    "    files_fov_id = df_fov_file_assign.loc[df_fov_file_assign['fov_id'] == fov_id]\n",
    "\n",
    "    \n",
    "    print('  spor : ')\n",
    "\n",
    "    files_spor = files_fov_id.loc[files_fov_id['cell_type'] == 'spor']\n",
    "    \n",
    "    \n",
    "    \n",
    "    print('    {} : '.format(protein1), end = '')\n",
    "    \n",
    "    protein1_spor_name = files_spor.loc[files_spor['channel'] == protein1]['filename']\n",
    "\n",
    "    if protein1_spor_name.empty:\n",
    "        print('--')\n",
    "    else:\n",
    "        print(protein1_spor_name.iloc[0])\n",
    "      \n",
    "    \n",
    "    \n",
    "    print('    {} : '.format(protein2), end = '')\n",
    "   \n",
    "    protein2_spor_name = files_spor.loc[files_spor['channel'] == protein2]['filename']\n",
    "\n",
    "    if protein2_spor_name.empty:\n",
    "        print('--')\n",
    "    else:\n",
    "        print(protein2_spor_name.iloc[0])\n",
    "    \n",
    "    \n",
    "        \n",
    "    print('  veg : ')\n",
    "\n",
    "    files_veg = files_fov_id.loc[files_fov_id['cell_type'] == 'veg']\n",
    "    \n",
    "    print('    {} : '.format(protein1), end = '')\n",
    "        \n",
    "    protein1_veg_name = files_veg.loc[files_veg['channel'] == protein1]['filename']\n",
    "\n",
    "    if protein1_veg_name.empty:\n",
    "        print('--')\n",
    "    else:\n",
    "        print(protein1_veg_name.iloc[0])\n",
    "        \n",
    "        \n",
    "            \n",
    "    print('    {} : '.format(protein2), end = '')\n",
    "\n",
    "    protein2_veg_name = files_veg.loc[files_veg['channel'] == protein2]['filename']\n",
    "\n",
    "    if protein2_veg_name.empty:\n",
    "        print('--')\n",
    "    else:\n",
    "        print(protein2_veg_name.iloc[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Main loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "129fd0ec8aaa4547a85f52294ab87323",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing FOVs:   0%|          | 0/1 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No file for channel 1 for spor cells in FOV 1.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1eaf26ff102740a499da072128f5c97d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing picks:   0%|          | 0/8 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculation finished.\n",
      "Total: 8 rings analysed, 0 excluded.\n"
     ]
    }
   ],
   "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",
    "for fov_id in tqdm(range(1, int(df_fov_file_assign['fov_id'].max())+1), desc=\"Processing FOVs\"):\n",
    "    for cell_type in ['spor', 'veg']:\n",
    "        \n",
    "        filename1 = df_fov_file_assign.loc[(df_fov_file_assign['fov_id']==fov_id)\n",
    "                                              &(df_fov_file_assign['cell_type']==cell_type)\n",
    "                                              &(df_fov_file_assign['channel']==protein1)]['filename']\n",
    "        filename2 = df_fov_file_assign.loc[(df_fov_file_assign['fov_id']==fov_id)\n",
    "                                              &(df_fov_file_assign['cell_type']==cell_type)\n",
    "                                              &(df_fov_file_assign['channel']==protein2)]['filename']\n",
    "        \n",
    "        if filename1.empty:\n",
    "            print('No file for channel 1 for {} cells in FOV {}.'.format(cell_type, str(fov_id)))\n",
    "            continue\n",
    "        if filename2.empty:\n",
    "            print('No file for channel 2 for {} cells in FOV {}.'.format(cell_type, str(fov_id)))\n",
    "            continue\n",
    "    \n",
    "        filename1 = filename1.iloc[0]\n",
    "        filename2 = filename2.iloc[0]\n",
    "\n",
    "        #load locs and convert distances from px to nm (Attention!)\n",
    "        locs1, info1, pixelsize = load_data(os.path.join(path,filename1))        \n",
    "        locs2, info2, pixelsize = load_data(os.path.join(path,filename2))\n",
    "        \n",
    "        #print(filename1)\n",
    "        #print(filename2)\n",
    "        # iterate over picks in a file\n",
    "        for pick in tqdm(np.unique(locs1.group), desc=\"Processing picks\"):\n",
    "\n",
    "            # select locs from pick\n",
    "            pick_locs1 = locs1[locs1.group == pick]\n",
    "            pick_locs2 = locs2[locs2.group == pick]\n",
    "        \n",
    "        \n",
    "            ################\n",
    "            # (1) estimate the postition of the ring using a histogram along the pick direction\n",
    "            #  a ring should yield a gaussian peak (2D projection in XY)\n",
    "            ################\n",
    "\n",
    "            #print('pick ID', pick)\n",
    "            #print('channel 1', len(pick_locs1))\n",
    "            #print('channel 2', len(pick_locs2))\n",
    "\n",
    "            max_bin = max(pick_locs1.y_pick_rot.max(), pick_locs2.y_pick_rot.max())\n",
    "            step_size = max_bin/binning\n",
    "            bins = np.arange(0,max_bin+1, step_size)\n",
    "            estimated_peak1, r_par1, hist_data1 = find_peak(pick_locs1, binning=binning, bins = bins, axes=\"y\")\n",
    "            estimated_peak2, r_par2, hist_data2 = find_peak(pick_locs2, binning=binning, bins = bins, 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",
    "\n",
    "            # seperate ring parameter to the corresponding ring\n",
    "            ring_parameter1 = [] # to work as the double band script\n",
    "            ring_parameter1.append(r_par1)\n",
    "            ring_parameter2 = [] # to work as the double band script\n",
    "            ring_parameter2.append(r_par2)\n",
    "\n",
    "            \n",
    "            peak_dist = ring_parameter1[0][1]-ring_parameter2[0][1]\n",
    "            \n",
    "            \n",
    "            \n",
    "            # set up plot with gridspec\n",
    "            fig = plt.figure(figsize=(18, 14), constrained_layout=True)\n",
    "            gs = fig.add_gridspec(2,2)\n",
    "            fig.suptitle((\"FOV {}, {}, Pick {} - Ring analysis\\n\"\n",
    "                          \"{}: {}\\n\"\n",
    "                          \"{}: {}\").format(fov_id, cell_type, pick, protein1, filename1, protein2, filename2),\n",
    "                          fontsize=16,\n",
    "                          ha=\"center\")\n",
    "            \n",
    "            ax1 = fig.add_subplot(gs[0, 0])\n",
    "           \n",
    "            \n",
    "            plot_peak_dist(pick_locs1,\n",
    "                           hist_data1,\n",
    "                           estimated_peak1,\n",
    "                           r_par1,\n",
    "                           binning=binning,\n",
    "                           cutoff=cutoff,\n",
    "                           axes=\"y_pick_rot\",\n",
    "                           title=protein1,\n",
    "                           color_hist = blue, \n",
    "                           color_fit = blue,\n",
    "                           ax=ax1)\n",
    "            \n",
    "            ax2 = fig.add_subplot(gs[1, 0])\n",
    "\n",
    "            plot_peak_dist(pick_locs2,\n",
    "                           hist_data2,\n",
    "                           estimated_peak2,\n",
    "                           r_par2,\n",
    "                           binning=binning,\n",
    "                           cutoff=cutoff,\n",
    "                           axes=\"y_pick_rot\",\n",
    "                           title=protein2,\n",
    "                           color_hist = orange, \n",
    "                           color_fit = orange,\n",
    "                           ax=ax2)\n",
    "            \n",
    "            ax3 = fig.add_subplot(gs[1,1])\n",
    "            \n",
    "            plot_peak_dist(pick_locs1,\n",
    "                           hist_data1,\n",
    "                           estimated_peak1,\n",
    "                           r_par1,\n",
    "                           binning=binning,\n",
    "                           cutoff=cutoff,\n",
    "                           axes=\"y_pick_rot\",\n",
    "                           title='Distance from {} to {}: {:.1f} nm'.format(protein1, protein2, peak_dist),\n",
    "                           color_hist = blue, \n",
    "                           color_fit = blue,\n",
    "                           annotate = False,\n",
    "                           ax=ax3)\n",
    "            \n",
    "            plot_peak_dist(pick_locs2,\n",
    "                           hist_data2,\n",
    "                           estimated_peak2,\n",
    "                           r_par2,\n",
    "                           binning=binning,\n",
    "                           cutoff=cutoff,\n",
    "                           axes=\"y_pick_rot\",\n",
    "                           title='Distance from {} to {}: {:.1f} nm'.format(protein1, protein2, peak_dist),\n",
    "                           color_hist = orange, \n",
    "                           color_fit = orange,\n",
    "                           annotate = False,\n",
    "                           ax=ax3)\n",
    "            \n",
    "            \n",
    "            # plot 3D ring\n",
    "            ax4 = fig.add_subplot(gs[0, 1], projection=\"3d\")\n",
    "            plot_3d_ring(pick_locs1, pick_locs2, blue, orange, protein1, protein2, ax=ax4)\n",
    "\n",
    "            \n",
    "            ax2.set_xlim(ax1.get_xlim())\n",
    "            ax3.set_xlim(ax1.get_xlim())\n",
    "\n",
    "                    \n",
    "            \n",
    "            ################\n",
    "            # (5) append data into large array\n",
    "            ################\n",
    "            \n",
    "            ring_data.append([fov_id, #running file index\n",
    "                              cell_type, # veg or spore\n",
    "                              protein1,\n",
    "                              protein2,\n",
    "                              filename1, #filename of protein species1\n",
    "                              filename2, #filename of protein species2\n",
    "                              pick, #pick number\n",
    "                              peak_dist, #peak-to-peak distance (nm)\n",
    "                              ])\n",
    "            \n",
    "            ################\n",
    "            # (6) save plots\n",
    "            ################\n",
    "            \n",
    "            # ring analysis\n",
    "            img_fname = \"fov_{}_{}_pick_{}\".format(fov_id, cell_type, pick)\n",
    "            img_name = os.path.join(analysis_folder, img_fname)\n",
    "            fig.savefig(img_name+img_format, dpi=dpi, format=\"png\")\n",
    "            plt.close(fig)\n",
    "            \n",
    "            \n",
    "       \n",
    "        \n",
    "print(\"Calculation finished.\")\n",
    "print(\"Total: {} rings analysed, {} excluded.\".format(len(ring_data), len(rings_excluded)))\n",
    "\n",
    "\n",
    "\n",
    "# Save analysis data\n",
    "df_ring_data_add = pd.DataFrame(ring_data, columns=[\"fov_id\",\n",
    "                                                \"cell_type\",\n",
    "                                                \"protein1\",\n",
    "                                                \"protein1\",\n",
    "                                                \"filename1\",\n",
    "                                                \"filename2\",\n",
    "                                                \"group\",\n",
    "                                                \"p2p distance\"\n",
    "                                                ])\n",
    "if df_ring_data is not None:\n",
    "    df_ring_data = pd.concat([df_ring_data, df_ring_data_add], ignore_index = True)\n",
    "else:\n",
    "    df_ring_data = df_ring_data_add\n",
    "    \n",
    "df_ring_data.to_csv(\"exchange_data.csv\")\n",
    "# Save dataframe with cell means for easy loading of data for postprocessing\n",
    "df_ring_data.to_pickle(\"exchange_data.pkl\")\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Results of previously analyzed datasets were detected.\n"
     ]
    }
   ],
   "source": [
    "df_ring_data = load_ring_data_df(path, \"exchange_data.pkl\")\n",
    "\n",
    "data_veg = df_ring_data[df_ring_data['cell_type'] == 'veg']\n",
    "data_spor = df_ring_data[df_ring_data['cell_type'] == 'spor']\n",
    "\n",
    "# histogram\n",
    "max_bin = max(data_veg['p2p distance'].max(), data_spor['p2p distance'].max())\n",
    "min_bin = min(data_veg['p2p distance'].min(), data_spor['p2p distance'].min())\n",
    "\n",
    "N_bins = 20\n",
    "step_size = (max_bin-min_bin)/N_bins\n",
    "bins = np.arange(min_bin,max_bin+1, step_size)\n",
    "\n",
    "n_veg, bins_veg = np.histogram(data_veg['p2p distance'], bins=bins)\n",
    "centers_veg = (bins[:-1] + bins[1:]) / 2\n",
    "\n",
    "n_spor, bins_spor = np.histogram(data_spor['p2p distance'], bins=bins)\n",
    "centers_spor = (bins[:-1] + bins[1:]) / 2\n",
    "\n",
    "binwidth = bins[1]-bins[0]\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot()\n",
    "\n",
    "ax.bar(centers_spor, n_spor, width=binwidth, color=red, alpha = 0.5, label = 'spor')\n",
    "ax.bar(centers_veg, n_veg, width=binwidth, color=blue, alpha = 0.5, label = 'veg')\n",
    "\n",
    "ax.axvline(0,c='gray', linewidth=2, linestyle=\"--\")\n",
    "\n",
    "ax.set_title('Distance from {} to {}'.format(protein1, protein2),fontsize=14)\n",
    "ax.set_xlabel(\"Distance (nm)\")\n",
    "ax.set_ylabel(\"Counts\")  \n",
    "    \n",
    "ax.legend()\n",
    "\n",
    "# ring analysis\n",
    "img_fname = \"Distance\"\n",
    "img_name = os.path.join(analysis_folder, img_fname)\n",
    "fig.savefig(img_name+img_format, dpi=dpi, format=\"png\")\n",
    "plt.close(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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
}