{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "by Andreas Bauer and David Boehringer, 19.06.2020\n",
    "This script contains two examples for 3D-PIV: The shift of a bar of binary pixels in one direction, the expansion\n",
    "and a real data set where we recorded two stacks of collagen fibres at the same field of view with confocal microscopy\n",
    " in reflection mode.  One stack contains a NK cell that deforms the matrix and the other doe not.\n",
    "Please download the data at https://github.com/fabrylab/3D_piv_example_data.git (180 MB, unpacked) and provide the\n",
    "folder in the code below.\n",
    "We tested this on ubuntu 16 and 18, with Anaconda Python installation. The whole script\n",
    "takes about 5 minutes on my 4 core-intel i5 @2.5 GHz Laptop. You should have !!! 8 Gb ob Memory !!!! or take care not\n",
    "to open  all matplotlib plots as interactive windows at once.\n",
    "For questions contact andreas.b.bauer@fau.de\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openpiv.pyprocess3D import *\n",
    "from openpiv.PIV_3D_plotting import *\n",
    "from openpiv.validation import sig2noise_val\n",
    "from openpiv.filters import replace_outliers\n",
    "from openpiv.lib import replace_nans\n",
    "import glob as glob\n",
    "import os\n",
    "from natsort import natsorted\n",
    "import matplotlib.animation as animation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make save_plots = True if you want to compare the \n",
    "visual results "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "save_plots = False\n",
    "out_put_folder = \"output_3D_test\"\n",
    "if save_plots:\n",
    "    if not os.path.exists(out_put_folder):\n",
    "        try:\n",
    "            os.mkdir(out_put_folder)\n",
    "        except:\n",
    "            print(\"could not generate output folder\")\n",
    "            save_plots = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "############ a group of bars shifted by 1 pixel to the each dimesion the second frame #############\n",
    "takes ~4 seconds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# constructing frame 1 and frame 2\n",
    "size = (32, 32, 32)\n",
    "shape1 = np.zeros(size)\n",
    "shape2 = np.zeros(size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shape1[16, 16, 25:27] = 1\n",
    "shape1[16, 16, 7:9] = 1\n",
    "shape1[16, 25:27, 16] = 1\n",
    "shape1[16, 7:9, 16] = 1\n",
    "shape1[25:27, 16, 16] = 1\n",
    "shape1[7:9, 16, 16] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "shape2[16, 16, 24:26] = 1\n",
    "shape2[16, 16, 8:10] = 1\n",
    "shape2[16, 24:26, 16] = 1\n",
    "shape2[16, 8:10, 16] = 1\n",
    "shape2[24:26, 16, 16] = 1\n",
    "shape2[8:10, 16, 16] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "window_size = (4, 4, 4)\n",
    "overlap = (3, 3, 3)\n",
    "search_area = (5, 5, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u, v, w, sig2noise = extended_search_area_piv3D(shape1, shape2, window_size=window_size, overlap=overlap,\n",
    "                                                search_area_size=search_area, subpixel_method='gaussian',\n",
    "                                                sig2noise_method='peak2peak',\n",
    "                                                width=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pdb\n",
    "# displaying the shapes with 3D scatter plot\n",
    "fig1 = scatter_3D(shape1, control=\"size\")\n",
    "fig2 = scatter_3D(shape2, control=\"size\")\n",
    "# 3d plot of the signal-to-noise rations\n",
    "fig3 = scatter_3D(sig2noise, control=\"size\")\n",
    "# 3d quiver plot of the displacement field\n",
    "fig4 = quiver_3D(-u, v, w, cmap=\"coolwarm\", quiv_args={\"arrow_length_ratio\":0.6})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# saving the plots\n",
    "if save_plots:\n",
    "    fig1.savefig(os.path.join(out_put_folder, \"displaced_bar_frame1.png\"))\n",
    "    fig2.savefig(os.path.join(out_put_folder, \"displaced_bar_frame2.png\"))\n",
    "    fig3.savefig(os.path.join(out_put_folder, \"displaced_bar_sig2noise.png\"))\n",
    "    fig4.savefig(os.path.join(out_put_folder, \"displaced_bar_deformation_field.png\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "################### test to check the replace_nans_function ######################\n",
    "takes ~4 seconds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ball shape with a gap of nans in the middle\n",
    "center = (5, 5, 5)\n",
    "size = (10, 10, 10)\n",
    "distance = np.linalg.norm(np.subtract(np.indices(size).T, np.asarray(center)), axis=len(center))\n",
    "arr = np.ones(size) * (distance <= 5)\n",
    "hide = arr == 0\n",
    "arr[5:7] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# displaying in 3d plots. Values outside of the original ball are hidden by setting to nan\n",
    "arr_show = arr.copy()\n",
    "arr_show[hide] = np.nan\n",
    "fig9 = scatter_3D(arr_show, size=50, sca_args={\"alpha\": 0.6})\n",
    "# replacing outliers\n",
    "arr = replace_nans(arr, max_iter=2, tol=2, kernel_size=2, method='disk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# displaying in 3d plots. Values outside of the original ball are hidden by setting to nan\n",
    "arr_show = arr.copy()\n",
    "arr_show[hide] = np.nan\n",
    "fig10 = scatter_3D(arr_show, size=50, sca_args={\"alpha\": 0.6})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# saving the plots\n",
    "if save_plots:\n",
    "    fig9.savefig(os.path.join(out_put_folder, \"replace_nan_gap.png\"))\n",
    "    fig10.savefig(os.path.join(out_put_folder, \"replace_nan_filled.png\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#################### real data example ############################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we recorded stacks of collagen fibres with confocal microscopy in reflection mode\n",
    "\"alive\" stack contains a force generating NK-cell, marked by the red circle in the animation\n",
    "\"relaxed\" stack is the same field of view with out the cell\n",
    "download the data at https://github.com/fabrylab/3D_piv_example_data.git\n",
    "this calculation takes ~ 3-4 minutes on my 4-core Intel i5@2.5 GHz Laptop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# please enter the path to the dataset provided at\n",
    "folder = r\"test_3d\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(folder): \n",
    "    import git \n",
    "    repo = git.Repo.clone_from(\"https://github.com/fabrylab/3D_piv_example_data.git\", './test_3d', branch='master')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(folder):\n",
    "    raise FileNotFoundError(\"path to 3d piv data '%s' does not exists\\n\"\n",
    "                            \". Please download the data from https://github.com/fabrylab/3D_piv_example_data.git\" % folder)\n",
    "# stack properties\n",
    "# factors for voxel size\n",
    "du = 0.2407\n",
    "dv = 0.2407\n",
    "dw = 1.0071\n",
    "# total image dimension for x y z\n",
    "image_dim = (123.02, 123.02, 122.86)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# keep these values for our nk cells stacks\n",
    "win_um = 12  # window size in µm\n",
    "fac_overlap = 0.3  # overlap in percent of the window size\n",
    "signoise_filter = 1.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# window size for stacks in pixel\n",
    "window_size = (int(win_um / du), int(win_um / dv), int(win_um / dw))\n",
    "overlap = (int(fac_overlap * win_um / du), int(fac_overlap * win_um / dv), int(fac_overlap * win_um / dw))\n",
    "search_area = (int(win_um / du), int(win_um / dv), int(win_um / dw))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load tense stacks\n",
    "images = natsorted(glob.glob(os.path.join(folder, \"Series001_t22_z*_ch00.tif\")))\n",
    "im_shape = plt.imread(images[0]).shape\n",
    "alive = np.zeros((im_shape[0], im_shape[1], len(images)))\n",
    "for i, im in enumerate(images):\n",
    "    alive[:, :, i] = plt.imread(im)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load relaxed stack\n",
    "images = natsorted(glob.glob(os.path.join(folder, \"Series003_t05_z*_ch00.tif\")))\n",
    "im_shape = plt.imread(images[0]).shape\n",
    "relax = np.zeros((im_shape[0], im_shape[1], len(images)))\n",
    "for i, im in enumerate(images):\n",
    "    relax[:, :, i] = plt.imread(im)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3D PIV\n",
    "u, v, w, sig2noise = extended_search_area_piv3D(relax, alive, window_size=window_size, overlap=overlap,\n",
    "                                                search_area_size=search_area, dt=(1 / du, 1 / dv, 1 / dw),\n",
    "                                                subpixel_method='gaussian',\n",
    "                                                sig2noise_method='peak2peak',\n",
    "                                                width=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# correcting stage drift between the field of views\n",
    "u -= np.nanmean(u)\n",
    "v -= np.nanmean(v)\n",
    "w -= np.nanmean(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# filtering\n",
    "uf, vf, wf, mask = sig2noise_val(u, v, w=w, sig2noise=sig2noise, threshold=signoise_filter)\n",
    "uf, vf, wf = replace_outliers(uf, vf, wf, max_iter=1, tol=100, kernel_size=2, method='disk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plotting\n",
    "# representation of the image stacks by maximums projections. The red circle marks the position of the cell\n",
    "def update_plot(i, ims, ax):\n",
    "    a1 = ax.imshow(ims[i])\n",
    "    a2 = ax.add_patch(plt.Circle((330, 140), 100, color=\"red\", fill=False))\n",
    "    return [a1, a2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ims = [np.max(relax[:, :, 60:], axis=2), np.max(alive[:, :, 60:], axis=2)]\n",
    "fig = plt.figure()\n",
    "ax = plt.gca()\n",
    "ani = animation.FuncAnimation(fig, update_plot, 2, interval=200, blit=True, repeat_delay=0, fargs=(ims, ax))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# unfiltered 3d deformation field\n",
    "fig11 = quiver_3D(-u, v, w, quiv_args={\"length\": 2, \"alpha\": 0.8, \"linewidth\": 1}, filter_def=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# filtered 3d deformation field\n",
    "fig12 = quiver_3D(-uf, vf, wf, quiv_args={\"length\": 2, \"alpha\": 0.8, \"linewidth\": 1},\n",
    "                  filter_def=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# saving the plots\n",
    "if save_plots:\n",
    "    fig11.savefig(os.path.join(out_put_folder, \"real_data_unfiltered.png\"))\n",
    "    fig12.savefig(os.path.join(out_put_folder, \"real_data_filtered.png\"))\n",
    "\n",
    "    # This needs a working ImageMagick installation, and probably works only on linux\n",
    "    try:\n",
    "        import imageio\n",
    "        plt.ioff()\n",
    "        f1 = plt.figure()\n",
    "        plt.imshow(ims[0])\n",
    "        plt.gca().add_artist(plt.Circle((330, 140), 100, color=\"red\", fill=False))\n",
    "        f1.savefig(os.path.join(out_put_folder,\"tem1.png\"))\n",
    "\n",
    "        f2 = plt.figure()\n",
    "        plt.imshow(ims[1])\n",
    "        plt.gca().add_artist(plt.Circle((330, 140), 100, color=\"red\", fill=False))\n",
    "        f2.savefig(os.path.join(out_put_folder,\"tem2.png\"))\n",
    "\n",
    "        i1 = plt.imread(os.path.join(out_put_folder,\"tem1.png\"))\n",
    "        i2 = plt.imread(os.path.join(out_put_folder, \"tem2.png\"))\n",
    "        imageio.mimsave(os.path.join(out_put_folder, \"reaL_data_max_proj.gif\"),[i1,i2], fps=1)\n",
    "        os.remove(os.path.join(out_put_folder,\"tem1.png\"))\n",
    "        os.remove(os.path.join(out_put_folder,\"tem2.png\"))\n",
    "        plt.ion()\n",
    "    except Exception as e:\n",
    "        print (\"failed to write gif of collagen embedded cell:\")\n",
    "        print(e)"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "formats": "ipynb,py:percent"
  },
  "kernelspec": {
   "display_name": "Python [conda env:openpiv] *",
   "language": "python",
   "name": "conda-env-openpiv-py"
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}