Customized image segmentation algorithm to find a sphere using cross-correlation and template matching

custom_feature_detection_template_matching.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "bbe81409-b19e-42f4-8c96-65e028a8be58",
    "deepnote_cell_type": "code"
   },
   "source": "##### CHANGEABLE PARAMETERS ######",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00001-8f342c55-853f-4f08-8031-68dc7309bb7f",
    "deepnote_cell_type": "code"
   },
   "source": "import os # to know which computer it is\nimport yaml",
   "execution_count": 2,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00002-fd40fc37-77fb-42c3-b3c6-438dbfabd9ab",
    "deepnote_cell_type": "code"
   },
   "source": "# we do not need to change Jupyter notebook contents for every computer\n# better to ignore this file from git probably, or be ready to commit it\n# every time we switch the folder\n\nyaml_file = 'path_location.yaml'\nyaml_args = yaml.load(open(yaml_file,'r'),yaml.CLoader)\npath_args = yaml_args['path']\nframe_ranges = yaml_args['frame_ranges']\nxlim = yaml_args['xlim']\n\nfor key in yaml_args:\n    print(key, yaml_args[key])",
   "execution_count": 3,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": "path G:\\Shared drives\\Liron_Simon_Thesis\\Experimental_data\\Raw_data\\060222\\C001H001S0001\nframe_ranges [160, 560]\nxlim [90, 150]\nfps 250\npixel2meter 0.126/1024\n"
    }
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00003-09d42e0d-e3ad-4a23-8d5f-889a9ffddb3d",
    "deepnote_cell_type": "code"
   },
   "source": "# date_taken = '150222'\n# file = '0001'\ndate_taken = path_args.split(os.sep)[-2]\nfile = path_args.split(os.sep)[-1].split('S')[-1]\n\n\nsigma = 1 # gaussian\nN = 5 # moving average window for the background\n\n# data position and additional parameters\n\n# DATA_PATH = f\"/media/ron/Liron USB/06-02-22/C001H001S{file}/\"+\"*.tif\" #\"../../C001H001S0001/*.tif\"\n# if os.environ['COMPUTERNAME'] == 'DESKTOP-OV3L1FT':\n#     DATA_PATH = r\"C:\\Users\\alex\\Downloads\\C001H001S0002\\*.tif\"\n# elif os.environ['COMPUTERNAME'] == 'DESKTOP-S68P8PR': # OFFICE PC\n#     DATA_PATH = rf\"G:\\Shared drives\\Liron_Simon_Thesis\\Experimental data\\Raw data\\{date_taken}\\C001H001S{file}\\S{file}\"\n\n# see link to Google Drive https://docs.google.com/spreadsheets/d/1dkG03RMQC-XQx8klTu3SgG6a0SPWhHm2aTcwBeC51Zc/edit#gid=1532686687&range=B7\n\n\nDATA_PATH = os.path.join(path_args,'*.tif')\n\n",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00004-a700c8f7-1607-4008-949c-9a71a174485e",
    "deepnote_cell_type": "code"
   },
   "source": "import numpy as np\nimport pandas as pd\nimport trackpy as tp\nimport pims\n\nfrom glob import glob\nfrom skimage.io.collection import alphanumeric_key\n\nfrom skimage.feature import match_template\nfrom skimage.registration import phase_cross_correlation\nfrom skimage.feature import peak_local_max\nfrom skimage.filters import gaussian\n",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00005-b0965d32-693b-425f-8710-ca20147907b1",
    "deepnote_cell_type": "code"
   },
   "source": "import matplotlib.pyplot as plt\n#matplotlib notebook # if you need interactivity like zoom\n%matplotlib inline\n\nimport matplotlib\nmatplotlib.rcParams['figure.figsize'] = (12, 12)",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00006-55c237ea-8725-4b63-ab64-3b40add99ee2",
    "deepnote_cell_type": "code"
   },
   "source": "filenames = sorted( glob(DATA_PATH), key=alphanumeric_key )\nprint(f\"Found {len(filenames)} files in this path {DATA_PATH} \")",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00007-bfb32bb7-6442-4c34-a6b2-f2925f42983b",
    "deepnote_cell_type": "code"
   },
   "source": "if len(frame_ranges) < 1:\n    frame_ranges = [0, len(filenames)]",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00008-ae2ba499-c12f-481d-8dff-a6fd9b7ee863",
    "deepnote_cell_type": "code"
   },
   "source": "# follow the example\n# http://soft-matter.github.io/trackpy/v0.5.0/tutorial/custom-feature-detection.html\n@pims.pipeline\ndef gray(fn):\n    # img = imread(fn)\n    return np.round(fn[:,slice(*xlim)]/4096*255).astype(np.uint8) # note the small ROI\n\n\n\nrawframes = gray(pims.open(os.path.join(DATA_PATH)))\nrawframes = rawframes[slice(*frame_ranges)]",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00009-edf70dca-591e-4921-9e37-cb241da1cb7e",
    "deepnote_cell_type": "code"
   },
   "source": "def moving_average(a, n=N) :\n    ret = np.cumsum(a, dtype=float, axis=0)\n    ret[n:,:,:] = ret[n:,:,:] - ret[:-n,:,:]\n    return ret[n - 1:,:,:] / n",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00010-ac69f31e-3175-4434-8284-46a507d66933",
    "deepnote_cell_type": "code"
   },
   "source": "# moving average backgrounds (list of backgrounds)\nbkg = moving_average(rawframes,n=N);",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00011-e80d22da-1d11-459e-b6ba-612841628b22",
    "deepnote_cell_type": "code"
   },
   "source": "# note the clipping - bright spots change\n# along the image, so we have to clip it\n# we could possibly also do some \n# histogram equaliztion\n\nframes = []\nfor f, b in zip(rawframes[N-1:], bkg):\n    frames.append(gaussian(np.clip(f - b, 0, 80),sigma=sigma))",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00012-ed91e57a-a65f-4fe8-8cb7-f65b7434b5ec",
    "deepnote_cell_type": "code"
   },
   "source": "template = frames[20][135:160,15:40]\nplt.imshow(template, cmap = 'gray')",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00013-c8ec9319-50bf-4029-9b76-e38ca04afff3",
    "deepnote_cell_type": "code"
   },
   "source": "img = frames[20]\n\nc = match_template(img, template)\n\nx, y = np.unravel_index(np.argmax(c), c.shape)\ntemplate_width, template_height = template.shape\nrect = plt.Rectangle((y, x), template_height, template_width, \n                     color='r', fc='none')\nplt.figure(num=None, figsize=(8, 6), dpi=80)\nplt.gca().add_patch(rect)\nplt.imshow(img);\n\n# plt.ylim(350,400)\n",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00014-0ab0eaa8-858a-407e-80e6-681eee437e04",
    "deepnote_cell_type": "code"
   },
   "source": "img = frames[100]\n\nc = match_template(img, template)\n\nplt.figure(num=None, figsize=(8, 6), dpi=80)\nfor x, y in peak_local_max(c, threshold_abs=0.85, \n                           exclude_border = 5):\n    rect = plt.Rectangle((y, x), template_height, template_width,\n                          color='r', fc='none')\n    plt.plot(y+template_height/2,x+template_width/2,'bo')\n    plt.gca().add_patch(rect)\n\n    small_img = img[x:x+template_width,y:y+template_height]\n    shift, error, diffphase = phase_cross_correlation(small_img, template,\n                                                  upsample_factor=20)\n    plt.plot(y+template_height/2+shift[1],x+template_width/2+shift[0],'r+',markersize=10)\n\nplt.imshow(img);\n\n# plt.ylim(240,260)\n# plt.ylim(300,400)\n# plt.gca().invert_yaxis()",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00015-6128e3db-f246-423a-b1d0-cd91a5f45645",
    "deepnote_cell_type": "code"
   },
   "source": "collect_features = []\nfor num, img in enumerate(frames):\n    c = match_template(img, template)\n    for x, y in peak_local_max(c, threshold_abs=0.85, \n                           exclude_border = 5):\n        small_img = img[x:x+template_width,\\\n                        y:y+template_height]\n        shift, error, diffphase = phase_cross_correlation(small_img, \\\n            template, upsample_factor=20)\n\n        # note the notation of features and images is different:\n        collect_features.append(\n            pd.DataFrame([{'x': y+template_height/2+shift[1],\n                            'y': x+template_width/2+shift[0],\n                            'frame': num\n                            },]))\n        \nfeatures = pd.concat(collect_features)\n",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00016-c28b5dac-fc0b-4ad5-82c1-ee9235aa2868",
    "deepnote_cell_type": "code"
   },
   "source": "tp.annotate(features[features.frame == 175], frames[175]);",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00017-80bc2547-a9f8-4f32-bcd6-925d85834e76",
    "deepnote_cell_type": "code"
   },
   "source": "features.to_csv(f'features{date_taken}_{file}.csv')",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00018-f4d3c06f-48b5-490f-a743-cac7fcab93a2",
    "deepnote_cell_type": "code"
   },
   "source": "import trackpy.predict\n\n@trackpy.predict.predictor\ndef predict(t1, particle):\n    velocity = np.array((0, -1)) # See fakeframe()\n    return particle.pos + velocity * (t1 - particle.t)\n\n#tr = pandas.concat(trackpy.link_df_iter(f, 0.5, predictor=predict))\ntr = trackpy.link(features, 45, predictor=predict, memory=5)\n",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00019-3a725e6b-28de-4a0d-8b63-100b866838f8",
    "deepnote_cell_type": "code"
   },
   "source": "tp.plot_traj(tr,label=True)",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00020-883b19fb-3781-4e85-ae8c-ead15f8302e0",
    "deepnote_cell_type": "code"
   },
   "source": "# remove particles that do not move vertically\n\nlst_num = []\nfor p in set(tr.particle):\n    t = tr[tr.particle == p]\n    if t.y.max() - t.y.min() > 50:\n        lst_num.append(p)\n    \nprint(lst_num)",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00021-eb3e5985-bbe4-4fec-8de6-c15ed58ebb4b",
    "deepnote_cell_type": "code"
   },
   "source": "tr = tr[tr.particle.isin(lst_num)]\ntp.plot_traj(tr,label=True)",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00022-1c5afbb0-73a1-4690-8b72-304a660b066f",
    "deepnote_cell_type": "code"
   },
   "source": "tr.to_csv(f'tr_{date_taken}_{file}.csv')",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00023-049460ea-0275-4ec5-bcf1-cfef7cc6b1e2",
    "deepnote_cell_type": "code"
   },
   "source": "# lets sort data\n#tm = tr2.reset_index(drop=True)\ntm = tr.sort_values(['particle', 'frame'], ignore_index=True)",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00024-66457a57-4c4b-4100-8b51-3ebd3b29f5f0",
    "deepnote_cell_type": "code"
   },
   "source": "from scipy.signal import savgol_filter\ndef calc_V(V_col, fps=fps, pixel2meter=pixel2meter):\n    # delta = np.gradient(V_col)*Pixel2meter\n    delta = savgol_filter(V_col, window_length=11, polyorder=2, deriv=1)\n    delta *= pixel2meter\n    return delta * fps",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00025-48e14806-4573-4f48-911f-1972200cef62",
    "deepnote_cell_type": "code"
   },
   "source": "tm['sec'] = tm.loc[:,'frame'] /fps\ntm['Vx'] = calc_V(tm.loc[:,'x'])\ntm['Vy'] = calc_V(tm.loc[:,'y'])",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00026-fc1251d0-a079-49ad-9fbe-e3f235e9b9b2",
    "deepnote_cell_type": "code"
   },
   "source": "#tm.plot(x=\"sec\", y=[\"Vx\",\"Vy\"])\n\ntm_i = tm[tm['particle'] == 0]\ntm_i = tm_i[5:-5]\ntm_i.plot(x=\"sec\", y=[\"Vx\",\"Vy\"])",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "cell_id": "00027-c54cfd66-4f13-4ca1-8e07-159521860c75",
    "deepnote_cell_type": "code"
   },
   "source": "fig,ax = plt.subplots(2,1)\n\nfor p in set(tm.particle):\n    t = tm[tm.particle == p]\n    t1 = t.drop(t[(t.Vy < .1) | (t.Vy > 0.2) ].index)\n    y = t1.y.values\n    Vy = t1.Vy.values\n    Vx = t1.x.values\n    fr = t1.frame.values\n    ax[0].plot(y,Vy/Vy[0]+p*.05,label=p)\n    ind = np.argmin(np.abs(y-800))\n    ax[0].text(y[ind],Vy[ind]/Vy[0]+p*.05,str(p),fontsize=14)\n    # t1.plot(x='y',y='Vy',ax=ax)\n\n    ax[1].plot(y,Vx,label=p)\n\n\nax[0].set_xlim([20,810])\nax[0].set_xlabel('$y$',fontsize=14)\nax[0].set_ylabel('$V_y/V_y(t=0)$ (shifted for clarity)',fontsize=14)\n\nax[1].set_xlim([20,810])\nax[1].set_xlabel('$y$',fontsize=14)\nax[1].set_ylabel('$V_x$',fontsize=14)\n\n# plt.legend()",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": "<a style='text-decoration:none;line-height:16px;display:flex;color:#5B5B62;padding:10px;justify-content:end;' href='https://deepnote.com?utm_source=created-in-deepnote-cell&projectId=83d975df-8c0a-4829-9055-240c8ceb60e4' target=\"_blank\">\n<img alt='Created in deepnote.com' style='display:inline;max-height:16px;margin:0px;margin-right:7.5px;' src='data:image/svg+xml;base64,PD94bWwgdmVyc2lvbj0iMS4wIiBlbmNvZGluZz0iVVRGLTgiPz4KPHN2ZyB3aWR0aD0iODBweCIgaGVpZ2h0PSI4MHB4IiB2aWV3Qm94PSIwIDAgODAgODAiIHZlcnNpb249IjEuMSIgeG1sbnM9Imh0dHA6Ly93d3cudzMub3JnLzIwMDAvc3ZnIiB4bWxuczp4bGluaz0iaHR0cDovL3d3dy53My5vcmcvMTk5OS94bGluayI+CiAgICA8IS0tIEdlbmVyYXRvcjogU2tldGNoIDU0LjEgKDc2NDkwKSAtIGh0dHBzOi8vc2tldGNoYXBwLmNvbSAtLT4KICAgIDx0aXRsZT5Hcm91cCAzPC90aXRsZT4KICAgIDxkZXNjPkNyZWF0ZWQgd2l0aCBTa2V0Y2guPC9kZXNjPgogICAgPGcgaWQ9IkxhbmRpbmciIHN0cm9rZT0ibm9uZSIgc3Ryb2tlLXdpZHRoPSIxIiBmaWxsPSJub25lIiBmaWxsLXJ1bGU9ImV2ZW5vZGQiPgogICAgICAgIDxnIGlkPSJBcnRib2FyZCIgdHJhbnNmb3JtPSJ0cmFuc2xhdGUoLTEyMzUuMDAwMDAwLCAtNzkuMDAwMDAwKSI+CiAgICAgICAgICAgIDxnIGlkPSJHcm91cC0zIiB0cmFuc2Zvcm09InRyYW5zbGF0ZSgxMjM1LjAwMDAwMCwgNzkuMDAwMDAwKSI+CiAgICAgICAgICAgICAgICA8cG9seWdvbiBpZD0iUGF0aC0yMCIgZmlsbD0iIzAyNjVCNCIgcG9pbnRzPSIyLjM3NjIzNzYyIDgwIDM4LjA0NzY2NjcgODAgNTcuODIxNzgyMiA3My44MDU3NTkyIDU3LjgyMTc4MjIgMzIuNzU5MjczOSAzOS4xNDAyMjc4IDMxLjY4MzE2ODMiPjwvcG9seWdvbj4KICAgICAgICAgICAgICAgIDxwYXRoIGQ9Ik0zNS4wMDc3MTgsODAgQzQyLjkwNjIwMDcsNzYuNDU0OTM1OCA0Ny41NjQ5MTY3LDcxLjU0MjI2NzEgNDguOTgzODY2LDY1LjI2MTk5MzkgQzUxLjExMjI4OTksNTUuODQxNTg0MiA0MS42NzcxNzk1LDQ5LjIxMjIyODQgMjUuNjIzOTg0Niw0OS4yMTIyMjg0IEMyNS40ODQ5Mjg5LDQ5LjEyNjg0NDggMjkuODI2MTI5Niw0My4yODM4MjQ4IDM4LjY0NzU4NjksMzEuNjgzMTY4MyBMNzIuODcxMjg3MSwzMi41NTQ0MjUgTDY1LjI4MDk3Myw2Ny42NzYzNDIxIEw1MS4xMTIyODk5LDc3LjM3NjE0NCBMMzUuMDA3NzE4LDgwIFoiIGlkPSJQYXRoLTIyIiBmaWxsPSIjMDAyODY4Ij48L3BhdGg+CiAgICAgICAgICAgICAgICA8cGF0aCBkPSJNMCwzNy43MzA0NDA1IEwyNy4xMTQ1MzcsMC4yNTcxMTE0MzYgQzYyLjM3MTUxMjMsLTEuOTkwNzE3MDEgODAsMTAuNTAwMzkyNyA4MCwzNy43MzA0NDA1IEM4MCw2NC45NjA0ODgyIDY0Ljc3NjUwMzgsNzkuMDUwMzQxNCAzNC4zMjk1MTEzLDgwIEM0Ny4wNTUzNDg5LDc3LjU2NzA4MDggNTMuNDE4MjY3Nyw3MC4zMTM2MTAzIDUzLjQxODI2NzcsNTguMjM5NTg4NSBDNTMuNDE4MjY3Nyw0MC4xMjg1NTU3IDM2LjMwMzk1NDQsMzcuNzMwNDQwNSAyNS4yMjc0MTcsMzcuNzMwNDQwNSBDMTcuODQzMDU4NiwzNy43MzA0NDA1IDkuNDMzOTE5NjYsMzcuNzMwNDQwNSAwLDM3LjczMDQ0MDUgWiIgaWQ9IlBhdGgtMTkiIGZpbGw9IiMzNzkzRUYiPjwvcGF0aD4KICAgICAgICAgICAgPC9nPgogICAgICAgIDwvZz4KICAgIDwvZz4KPC9zdmc+' > </img>\nCreated in <span style='font-weight:600;margin-left:4px;'>Deepnote</span></a>",
   "metadata": {
    "tags": [],
    "created_in_deepnote_cell": true,
    "deepnote_cell_type": "markdown"
   }
  }
 ],
 "nbformat": 4,
 "nbformat_minor": 4,
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.7 ('liron_danny')",
   "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.7"
  },
  "vscode": {
   "interpreter": {
    "hash": "fa416831bb8d2bce8670c6668d9a9f308649692bc2c0170a27935576e90f9346"
   }
  },
  "deepnote_notebook_id": "68e5a673-cc4c-4023-8305-fd22dc421700",
  "deepnote": {},
  "deepnote_execution_queue": []
 }
}