PSF_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a1a26d4d-88d8-4965-bfa7-905b22cf2637",
   "metadata": {},
   "source": [
    "# **Average Point Spread Function (PSF) distillation**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14c62bd3-7b28-4c3f-8122-350631e5b4dc",
   "metadata": {},
   "source": [
    "#### This notebook makes an average PSF out of a z-stack with acquired fluorescence beads. It will also perform statistical analysis over individual beads and averaged PSF. Accepted formats are ***czi*** and ***tif***."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a53f6ae-32c3-4b85-9ad4-820f8dc13af0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "import numpy as np\n",
    "from skimage.io import imread, imsave\n",
    "from pyclesperanto_prototype import imshow\n",
    "import pyclesperanto_prototype as cle\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from aicsimageio import AICSImage\n",
    "import os\n",
    "import tifffile\n",
    "import stackview\n",
    "import xml.etree.ElementTree as ET\n",
    "from psf_functions import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb5b1311-4e0e-45bb-a4a6-6347d6489e47",
   "metadata": {},
   "source": [
    "The main analysis is done using [pyclesperanto-prototype ](https://pypi.org/project/pyclesperanto-prototype/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2de63396-8082-4756-acdf-88bf67a6f496",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipyfilechooser import FileChooser\n",
    "fc = FileChooser()\n",
    "display(fc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "394db810-8d8d-44c6-9e5e-a4b3666031d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_path = fc.selected_path + '\\\\' + fc.selected_filename\n",
    "print(\"File path: \", file_path)\n",
    "#path = os.path.abspath(\"\")\n",
    "#print(path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "690c5bc9-29ce-4aaa-acb6-172dc017e696",
   "metadata": {},
   "outputs": [],
   "source": [
    "head, tail = os.path.split(file_path)\n",
    "\n",
    "results_path = head + '\\\\results'  \n",
    "if not os.path.exists(results_path):\n",
    "    os.makedirs(results_path)\n",
    "\n",
    "if tail.endswith('tif') or tail.endswith('czi'):\n",
    "    file_name = tail[:-4]\n",
    "else:\n",
    "    print('Format not supported')\n",
    "\n",
    "path_save = results_path + '\\\\' + file_name\n",
    "print(path_save)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0957a684-a47d-4c71-8773-acaf15534032",
   "metadata": {},
   "outputs": [],
   "source": [
    "if file_path[-3:] == 'czi':\n",
    "    # Get an AICSImage object\n",
    "    img = AICSImage(file_path)\n",
    "    data = img.data\n",
    "    root = img.metadata\n",
    "    bead_image = cle.asarray(data[0][0])\n",
    "    print(bead_image.shape)\n",
    "    pixels, metadata_ome = metadata_read(file_path)\n",
    "    print('Voxel size (x,y,z) in nm: ', pixels)\n",
    "if file_path[-3:] == 'tif':\n",
    "    bead_image = imread(file_path)\n",
    "    print(bead_image.shape)\n",
    "    pixels, metadata_ome = metadata_read(file_path)\n",
    "    print('Voxel size (x,y,z) in nm: ', pixels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b76b4236-ca28-4666-83e4-751023aaaf41",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "stackview.imshow(cle.maximum_x_projection(bead_image), colorbar=True)\n",
    "stackview.imshow(cle.maximum_y_projection(bead_image), colorbar=True)\n",
    "stackview.imshow(cle.maximum_z_projection(bead_image), colorbar=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbb5617e-35ad-4a10-83c7-35d540f29cf8",
   "metadata": {},
   "outputs": [],
   "source": [
    "stackview.slice(bead_image, continuous_update=True, zoom_factor = 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58f6d225-35c0-48a2-ac59-e747131692d1",
   "metadata": {},
   "source": [
    "For determining an average PSF, technically we can crop out all the individual beads, align them and then average the images. Therefore, we segment the objects and determine their center of mass."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "389b7a12-36b9-4438-8704-9d8a9f18507d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Segment objects\n",
    "label_image = cle.voronoi_otsu_labeling(bead_image)\n",
    "stackview.imshow(label_image, labels=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46a30e46-62ae-47a5-ab83-f7c85c2b7bb9",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# determine center of mass for each object\n",
    "stats = cle.statistics_of_labelled_pixels(bead_image, label_image)\n",
    "\n",
    "df = pd.DataFrame(stats)\n",
    "df[[\"mass_center_x\", \"mass_center_y\", \"mass_center_z\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "049659e6-f05b-4ff3-ac47-3b6d8fe7c4bf",
   "metadata": {},
   "source": [
    "# PSF averaging"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4dcf5865-34bb-43b7-9eb4-1678cde40d95",
   "metadata": {},
   "source": [
    "Next, we will iterate over the beads and crop them out by translating them into a smaller PSF image. We calculate FWHM for each bead at the same time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7db12fa5-f3b9-4c12-8dd3-d20c535472c5",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# configure size of future PSF image\n",
    "psf_radius = 20\n",
    "size = psf_radius * 2 + 1\n",
    "\n",
    "# initialize PSF\n",
    "single_psf_image = cle.create([size, size, size])\n",
    "avg_psf_image = cle.create([size, size, size])\n",
    "\n",
    "# initialize FWHM\n",
    "FWHM_x = []\n",
    "FWHM_y = []\n",
    "FWHM_zx = []\n",
    "FWHM_zy = []\n",
    "\n",
    "num_psfs = len(df)\n",
    "indx = []\n",
    "for index, row in df.iterrows():\n",
    "    x = row[\"mass_center_x\"]\n",
    "    y = row[\"mass_center_y\"]\n",
    "    z = row[\"mass_center_z\"]\n",
    "    \n",
    "    print(\"Bead\", index, \"at position\", x, y, z)\n",
    "    \n",
    "    # move PSF in right position in a smaller image\n",
    "    cle.translate(bead_image, single_psf_image, \n",
    "                  translate_x= -x + psf_radius,\n",
    "                  translate_y= -y + psf_radius,\n",
    "                  translate_z= -z + psf_radius)\n",
    "\n",
    "    [fwhm_x, fwhm_y, fwhm_zx, fwhm_zy] = FWHM(size, pixels, single_psf_image)\n",
    "    print('fwhm x: ', fwhm_x, 'nm', ' fwhm y: ', fwhm_y, 'nm', ' fwhm zx: ', fwhm_zx, 'nm', ' fwhm zy: ', fwhm_zy, 'nm')\n",
    "    \n",
    "    # visualize\n",
    "    fig, axs = plt.subplots(1,3)\n",
    "    axs[0].set_title(\"Bead \" + str(index))\n",
    "    #plt.title(\"Bead \" + str(index))\n",
    "    stackview.imshow(cle.maximum_x_projection(single_psf_image), plot=axs[0])\n",
    "    stackview.imshow(cle.maximum_y_projection(single_psf_image), plot=axs[1])\n",
    "    stackview.imshow(cle.maximum_z_projection(single_psf_image), plot=axs[2])\n",
    "\n",
    "    #indx.append(index)\n",
    "        # average\n",
    "    avg_psf_image = avg_psf_image + single_psf_image / num_psfs\n",
    "    \n",
    "    FWHM_x = np.append(FWHM_x,fwhm_x)\n",
    "    FWHM_y = np.append(FWHM_y,fwhm_y)\n",
    "    FWHM_zx = np.append(FWHM_zx,fwhm_zx) \n",
    "    FWHM_zy = np.append(FWHM_zy,fwhm_zy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca08dfdc-5634-411f-9bdf-2fc418afec28",
   "metadata": {},
   "outputs": [],
   "source": [
    "# mean and standard deviation of FWHM value\n",
    "FWHM_x_mean = round(np.mean(FWHM_x))\n",
    "FWHM_x_std = round(np.std(FWHM_x))\n",
    "FWHM_y_mean =  round(np.mean(FWHM_y))\n",
    "FWHM_y_std = round(np.std(FWHM_y))\n",
    "FWHM_zx_mean = round(np.mean(FWHM_zx))\n",
    "FWHM_zx_std = round(np.std(FWHM_zx))\n",
    "FWHM_zy_mean = round(np.mean(FWHM_zy))\n",
    "FWHM_zy_std = round(np.std(FWHM_zy))\n",
    "\n",
    "print('FWHM in x: ', FWHM_x_mean, '+-',FWHM_x_std, 'nm')\n",
    "print('FWHM in y: ', FWHM_y_mean, '+-',FWHM_y_std, 'nm')\n",
    "print('FWHM in zx: ', FWHM_zx_mean, '+-',FWHM_zx_std, 'nm')\n",
    "print('FWHM in zy: ', FWHM_zy_mean, '+-',FWHM_zy_std, 'nm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ad96e5d-4663-480a-92d8-dae46f81e1a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(nrows=1,ncols=3,subplot_kw={'xticks': [], 'yticks': []})\n",
    "stackview.imshow(cle.maximum_x_projection(avg_psf_image), plot=axs[0])\n",
    "stackview.imshow(cle.maximum_y_projection(avg_psf_image), plot=axs[1])\n",
    "stackview.imshow(cle.maximum_z_projection(avg_psf_image), plot=axs[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d36ccdbe-bc55-48f1-823c-711d484cf8a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "avg_psf_image.min(), avg_psf_image.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0aee9212-3667-48df-8a7c-81e320c8aab7",
   "metadata": {},
   "outputs": [],
   "source": [
    "normalized_psf = avg_psf_image / np.sum(avg_psf_image)\n",
    "\n",
    "stackview.imshow(normalized_psf, colorbar=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4082c3d-5fd8-4b3a-b926-b1a2490f778c",
   "metadata": {},
   "outputs": [],
   "source": [
    "normalized_psf.min(), normalized_psf.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28619f12-a292-4655-bab9-559ebe6e63c8",
   "metadata": {},
   "source": [
    "**Saving projection tif image of averaged PSF**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b75ebd3b-fe3a-4a96-9e4e-cd657f7df1fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.savefig(path_save + \"_projections.tif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34ee5164-9665-4e0c-86a9-ba5f4f1fc192",
   "metadata": {},
   "source": [
    "**Saving average PSF stack in ome.tif file**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97bd2e7c-0170-465d-962d-7902c4eda194",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save the numpy array as OME TIFF with metadata\n",
    "tifffile.imwrite(path_save + '_avgPSF_stack.ome.tiff', normalized_psf, metadata=metadata_ome)\n",
    "#imsave(path_save + '_avgPSF_stack.tif', normalized_psf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "764c9919-47da-4c71-9614-1fbd0574aad4",
   "metadata": {},
   "source": [
    "### Analysis of PSF"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7ba8da6-bdce-464e-b567-a047ba221d1a",
   "metadata": {},
   "source": [
    "Data [conversion](https://github.com/clEsperanto/pyclesperanto_prototype/blob/master/demo/interoperability/numpy.ipynb) from pyclesperanto to numpy. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77350f7e-ff28-47c9-b63d-2ff8aaa766d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "z_projection = np.asarray(cle.maximum_z_projection(avg_psf_image))\n",
    "x_projection = np.asarray(cle.maximum_x_projection(avg_psf_image))\n",
    "y_projection = np.asarray(cle.maximum_y_projection(avg_psf_image))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88fb1967-9d6a-4755-aaff-fdbaf925fed4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Indexes of the maximum in image array:\n",
    "ind_x = np.unravel_index(np.argmax(x_projection, axis=None), x_projection.shape)\n",
    "ind_y = np.unravel_index(np.argmax(y_projection, axis=None), y_projection.shape)\n",
    "ind_z = np.unravel_index(np.argmax(z_projection, axis=None), z_projection.shape)\n",
    "\n",
    "print('x: ', ind_x)\n",
    "print('y: ', ind_y)\n",
    "print('z: ', ind_z)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eef350f3-b8b8-4b3a-bcfb-6f4ed1efeca1",
   "metadata": {},
   "source": [
    "#### FWHM for z-axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2934422-6faa-45f2-af6b-244a9f1549ee",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "units_xy, units_z, units_xynew, units_znew = units(size, pixels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f752580a-75c4-44db-a849-4c4e8dd8ae48",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normalized curve\n",
    "zx_resolution = x_projection[ind_x[0],:]/avg_psf_image.max()\n",
    "zy_resolution = y_projection[:,ind_y[1]]/avg_psf_image.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d7e0462-82de-4c7a-b7e6-111b3df6dda9",
   "metadata": {},
   "source": [
    "Cubic Spline Interpolation from [Python Numerical Methods](https://pythonnumericalmethods.studentorg.berkeley.edu/notebooks/chapter17.03-Cubic-Spline-Interpolation.html) for z-axis in x and y planes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82ae38e7-6be7-449f-b53d-752ac83bdb6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "f_zx = CubicSpline(units_z,zx_resolution, bc_type='natural')\n",
    "f_zy = CubicSpline(units_z,zy_resolution, bc_type='natural')\n",
    "units_znew = np.linspace(units_z[0], units_z[-1], 1000)\n",
    "fzx_new = f_zx(units_znew)\n",
    "fzy_new = f_zy(units_znew)\n",
    "\n",
    "indices_zx = np.where(fzx_new > 0.5)[0].tolist()\n",
    "fwhm_zx = round(units_znew[indices_zx[-1]]-units_znew[indices_zx[0]])\n",
    "indices_zy = np.where(fzy_new > 0.5)[0].tolist()\n",
    "fwhm_zy = round(units_znew[indices_zy[-1]]-units_znew[indices_zy[0]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5685b25c-77ee-4997-9434-fbe14778cb7c",
   "metadata": {},
   "source": [
    "[matplotlib](https://matplotlib.org/stable/users/explain/quick_start.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "463f6bd1-7a1b-4325-84ca-a24eb0797aec",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "ax.plot(units_znew,fzx_new,'black')\n",
    "ax.plot(units_znew, fzy_new, 'black')\n",
    "ax.plot(units_z,zx_resolution, 'ro', label='zx')\n",
    "ax.plot(units_z,zy_resolution, 'bo', label='zy')\n",
    "ax.set_title('FWHM')\n",
    "ax.legend()\n",
    "plt.show()\n",
    "print('zx fwhm: ', fwhm_zx, 'nm')\n",
    "print('zy fwhm: ', fwhm_zy, 'nm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44af2e74-8df6-4263-a06c-98a678d8d0dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.savefig(path_save + \"_z_fwhm.svg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "486a73aa-acc6-43f7-b03e-ddad35a78ee4",
   "metadata": {},
   "source": [
    "#### FWHM for xy plane"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1f10912-fbf5-485f-ad68-09bf83ec1d61",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_resolution = z_projection[ind_z[0],:]/avg_psf_image.max()\n",
    "y_resolution = z_projection[:,ind_z[1]]/avg_psf_image.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e6ee3dc-9f5f-4511-816e-ec932a5308b9",
   "metadata": {},
   "source": [
    "Cubic Spline Interpolation from [Python Numerical Methods](https://pythonnumericalmethods.studentorg.berkeley.edu/notebooks/chapter17.03-Cubic-Spline-Interpolation.html) for xy (lateral) resolution (in x and y direction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a583ecb-1afa-4186-b849-82078bd2b994",
   "metadata": {},
   "outputs": [],
   "source": [
    "f_x = CubicSpline(units_xy,x_resolution, bc_type='natural')\n",
    "f_y = CubicSpline(units_xy,y_resolution, bc_type='natural')\n",
    "units_xynew = np.linspace(units_xy[0], units_xy[-1], 1000)\n",
    "fx_new = f_x(units_xynew)\n",
    "fy_new = f_y(units_xynew)\n",
    "\n",
    "indices_x = np.where(fx_new > 0.5)[0].tolist()\n",
    "fwhm_x = round(units_xynew[indices_x[-1]]-units_xynew[indices_x[0]])\n",
    "indices_y = np.where(fy_new > 0.5)[0].tolist()\n",
    "fwhm_y = round(units_xynew[indices_y[-1]]-units_xynew[indices_y[0]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49d09363-c310-4651-91dd-3dfcbd801bc0",
   "metadata": {},
   "source": [
    "[matplotlib](https://matplotlib.org/stable/users/explain/quick_start.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ce03abb-7ac7-482d-b81c-3aa2fad13cc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "ax.plot(units_xynew,fx_new,'black')\n",
    "ax.plot(units_xy,x_resolution, 'ro', label='x')\n",
    "ax.set_title('FWHM')\n",
    "ax.legend()\n",
    "plt.show()\n",
    "print('x fwhm: ', fwhm_x, 'nm') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fd6602b-11f4-472b-bdeb-f5c477f81497",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.savefig(path_save + \"_x_fwhm.svg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edae9772-ec69-444f-bf7f-260040dde444",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "ax.plot(units_xynew, fy_new, 'black')\n",
    "ax.plot(units_xy,y_resolution, 'bo', label='y')\n",
    "ax.set_title('FWHM')\n",
    "ax.legend()\n",
    "plt.show()\n",
    "print('y fwhm: ', fwhm_y, 'nm') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e16d1e4d-0d5a-47cc-b3c9-ee46ad277c73",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.savefig(path_save + \"_y_fwhm.svg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af5f5fe3-33b2-478a-849f-969893d9d9dd",
   "metadata": {},
   "source": [
    "#### Lateral and axial resolution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82e576df-449c-4358-9e75-c708275baf07",
   "metadata": {},
   "source": [
    "Choose which direction you want to use: \n",
    "x vs. y in lateral\n",
    "zx vs. zy in axial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "010cb0b3-8ff9-456e-9a42-8f8c2e03584f",
   "metadata": {},
   "outputs": [],
   "source": [
    "lateral = x_resolution\n",
    "axial = zy_resolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9bdabfa7-99c2-46a1-88d9-e21dc5e6f4b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "if np.all(lateral == x_resolution):\n",
    "    lateral_fit = fx_new\n",
    "    lateral_fwhm = fwhm_x\n",
    "    which_xy = '(x)'\n",
    "else:\n",
    "    lateral_fit = fy_new\n",
    "    lateral_fwhm = fwhm_y\n",
    "    which_xy = '(y)'\n",
    "if np.all(axial == zx_resolution):\n",
    "    axial_fit = fzx_new\n",
    "    axial_fwhm = fwhm_zx\n",
    "    which_z = '(zx)'\n",
    "else: \n",
    "    axial_fit = fzy_new\n",
    "    axial_fwhm = fwhm_zy\n",
    "    which_z = '(zy)'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df6a0f7a-90ff-43df-836a-6f2b46338981",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "ax.plot(units_xynew, lateral_fit,'black')\n",
    "ax.plot(units_xy, lateral, 'ro', label='lateral')\n",
    "ax.plot(units_znew, axial_fit, 'black')\n",
    "ax.plot(units_z, axial, 'bo', label='axial')\n",
    "ax.set_title('FWHM')\n",
    "ax.set_xlabel('Distance (nm)')\n",
    "ax.set_ylabel('Normalized intensity (a.u.)')\n",
    "ax.legend()\n",
    "plt.show()\n",
    "print('lateral ' + which_xy + ': ', lateral_fwhm, 'nm')\n",
    "print('axial ' + which_z + ': ', axial_fwhm, 'nm')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a82df5e-e750-45e4-93a0-a3a3f4fc764f",
   "metadata": {},
   "source": [
    "Values for lateral and axial resolution come from the averaged psf calculated above and these values has no error by default. \n",
    "We can use the FWHM values from each bead and calculate standard deviation (calculated above as well) for the mean FWHM value with an error. Nevertheless, those two value slightly differs."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec6f20db-a7ba-4a9a-b7bf-c4539f18f5a4",
   "metadata": {},
   "source": [
    "**Saving the graph in the file**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c798576c-0464-4d46-afe3-a594cafd6799",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.savefig(path_save + \"_fwhm.svg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e968eef1-f2d2-4144-ae8c-7a3f09f3f28e",
   "metadata": {},
   "source": [
    "**Saving the data in the csv file**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bedd3b0c-c473-40a3-8dd9-a2895c8c7217",
   "metadata": {},
   "outputs": [],
   "source": [
    "axis = [\"lateral\", \"axial\"]\n",
    "datatype = [\"data\", \"Cubic Spline\"]\n",
    "datapoints = {\"distance_lateral\": units_xy, \"lateral\": x_resolution,\n",
    "              \"distance_axial\": units_z, \"axial\": zy_resolution}            \n",
    "spline = {\"distance_lateral\": units_xynew, \"CS_lateral\": fx_new,\n",
    "                 \"distance_axial\": units_znew, \"CS_axial\": fzy_new}\n",
    "\n",
    "df_data = pd.DataFrame(datapoints)\n",
    "df_spline = pd.DataFrame(spline)\n",
    "\n",
    "df_data.to_csv(path_save + '_datapoints.csv')\n",
    "df_spline.to_csv(path_save + '_spline.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1827f7a9-4b86-41f7-8510-80518d13a721",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(path_save + '_fwhm.txt', 'w') as f:\n",
    "    f.write('FWHMs of ' + '%d' %len(df) + ' beads:' + '\\n' +\n",
    "            'FWHM in x: ' + '%d' %FWHM_x_mean + ' +- ' + '%d' %FWHM_x_std + ' nm' + '\\n' +\n",
    "            'FWHM in y: ' + '%d' %FWHM_y_mean + ' +- ' + '%d' %FWHM_y_std + ' nm' + '\\n' +    \n",
    "            'FWHM in zx: ' + '%d' %FWHM_zx_mean + ' +- ' + '%d' %FWHM_zx_std + ' nm' + '\\n' + \n",
    "            'FWHM in zy: ' + '%d' %FWHM_zy_mean + ' +- ' + '%d' %FWHM_zy_std + ' nm' + '\\n' + '\\n' +\n",
    "            'Averaged PSF out of ' + '%d' %len(df) + ' beads:' +  '\\n' + \n",
    "            'Lateral ' + which_xy + ': ' + '%d' %lateral_fwhm + ' nm' + '\\n' +\n",
    "            'Axial ' + which_z + ': ' + '%d' %axial_fwhm + ' nm' + '\\n') "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8486e42c-91af-49ae-b7d9-ca8cd6456d20",
   "metadata": {},
   "source": [
    "**Saving individual projection in tif files**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ce7b583-9966-447d-afd1-cd1dee3ba32f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "plt.imshow(cle.maximum_z_projection(normalized_psf),cmap='gray')\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "fig.savefig(path_save + \"_z_projection.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68e2c7fc-d1f3-4f7a-90b6-ac0a5f268fdc",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "plt.imshow(cle.maximum_y_projection(normalized_psf),cmap='gray')\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "fig.savefig(path_save + \"_y_projection.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bdbd1c7-92f8-41d1-8473-13401c936fac",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "plt.imshow(cle.maximum_x_projection(normalized_psf),cmap='gray')\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "fig.savefig(path_save + \"_x_projection.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42e1e71c-8b18-400c-b189-0518b5afa3db",
   "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.9.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}