OLD_extract_features.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: Qt5Agg\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Anaconda3\\envs\\wradlib0_11\\lib\\site-packages\\h5py\\__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n"
     ]
    }
   ],
   "source": [
    "%matplotlib\n",
    "import wradlib\n",
    "#import cv2\n",
    "import numpy as np\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "#from matplotlib import animation\n",
    "import matplotlib.patches as mpatches\n",
    "from matplotlib.collections import PatchCollection\n",
    "from scipy.ndimage import zoom\n",
    "import datetime\n",
    "import warnings\n",
    "warnings.simplefilter('once', DeprecationWarning)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import ndimage as ndi\n",
    "from skimage import feature"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read sample data\n",
    "\n",
    "Data is from the German Weather Service: the so called RY product represents rainfall intensity composite for the whole of Germany in 5 minute intervals. \n",
    "\n",
    "Spatial resolution: `1 x 1 km`; spatial extent: `900 x 900 km`.\n",
    "\n",
    "**Information required from user**\n",
    "\n",
    "- specify the directory `datadir` where you store the RY data (unpack the ry archives there).\n",
    "- select a specific interval by commenting/uncommenting the `dtimes` lines.\n",
    "- decide whether you need to reduce the resolution (downsize the image by a `downsizeby`) in order to avoid memory problems (this becomes relevant once you solve the 2D-adveciton equation...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set data directory\n",
    "datadir = \"data/ry\"\n",
    "\n",
    "# Original grid dimensions\n",
    "nx = 900\n",
    "ny = 900\n",
    "\n",
    "# pixel size (in meters)\n",
    "dx = 1000.\n",
    "dy = 1000.\n",
    "\n",
    "# Downsize by factor \"downsizeby\"\n",
    "#    downsizeby = 1 will leave the dimensions unchanged,\n",
    "#    but for a 900x900 km grid, downsizing might be \n",
    "#    required in order to avoid MemoryError\n",
    "downsizeby = 1\n",
    "\n",
    "# interval between observations (in seconds)\n",
    "interval = 300\n",
    "\n",
    "# Set time window\n",
    "##dtimes = wradlib.util.from_to(\"2008-06-02 17:00:00\", \"2008-06-02 19:00:00\", interval)\n",
    "##dtimes = wradlib.util.from_to(\"2015-04-26 17:00:00\", \"2015-04-26 19:00:00\", interval)\n",
    "##dtimes = wradlib.util.from_to(\"2015-03-29 17:00:00\", \"2015-03-29 19:00:00\", interval)\n",
    "dtimes = wradlib.util.from_to(\"2016-05-29 16:00:00\", \"2016-05-29 19:00:00\", interval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\src\\git\\heistermann\\wradlib\\wradlib\\trafo.py:127: RuntimeWarning: divide by zero encountered in log10\n",
      "  return 10. * np.log10(x)\n"
     ]
    }
   ],
   "source": [
    "# Compute grid dimensions and grid coordinates after resampling\n",
    "dx2, dy2 = dx*downsizeby, dy*downsizeby\n",
    "nx2, ny2 = int(nx/downsizeby), int(ny/downsizeby)\n",
    "\n",
    "X2, Y2 = np.meshgrid( np.arange(0,nx2*dx2, dx2), np.arange(0,ny2*dy2, dy2) )\n",
    "\n",
    "# Define container\n",
    "frames = np.zeros( (len(dtimes), nx2, ny2 ) )\n",
    "\n",
    "# Read the data, convert back to dBZ, and downsize\n",
    "#   (maybe also try with keeping mm/h instead of converting to dBZ?)\n",
    "for i, dtime in enumerate(dtimes):\n",
    "    fname = dtime.strftime( os.path.join(datadir, \"raa01-ry_10000-%y%m%d%H%M-dwd---bin\") )\n",
    "    frames[i] = zoom( wradlib.io.read_RADOLAN_composite(fname, missing=0)[0], 1./downsizeby, order=1)\n",
    "    frames[i] = wradlib.trafo.decibel( wradlib.zr.r2z(frames[i]) )\n",
    "    frames[i][frames[i]<0] = 0 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "#im = frames.sum(axis=0)\n",
    "im = frames[10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import measure\n",
    "from skimage import filters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Anaconda3\\envs\\wradlib0_11\\lib\\site-packages\\matplotlib\\cbook\\deprecation.py:106: MatplotlibDeprecationWarning: The spectral and spectral_r colormap was deprecated in version 2.0. Use nipy_spectral and nipy_spectral_r instead.\n",
      "  warnings.warn(message, mplDeprecation, stacklevel=1)\n"
     ]
    }
   ],
   "source": [
    "blobs = im > 5.\n",
    "\n",
    "all_labels = measure.label(blobs)\n",
    "blobs_labels = measure.label(blobs, background=0)\n",
    "\n",
    "plt.figure(figsize=(9, 3.5))\n",
    "plt.subplot(131)\n",
    "plt.imshow(blobs, cmap='gray')\n",
    "plt.axis('off')\n",
    "plt.subplot(132)\n",
    "plt.imshow(all_labels, cmap='spectral')\n",
    "plt.axis('off')\n",
    "plt.subplot(133)\n",
    "plt.imshow(blobs_labels, cmap='spectral')\n",
    "plt.axis('off')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "886"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nb_labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import ndimage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Anaconda3\\envs\\wradlib0_11\\lib\\site-packages\\scipy\\ndimage\\measurements.py:272: DeprecationWarning: In future, it will be an error for 'np.bool_' scalars to be interpreted as an index\n",
      "  return _nd_image.find_objects(input, max_label)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(-0.5, 34.5, 63.5, -0.5)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mask = frames[10] > 10\n",
    "\n",
    "label_im, nb_labels = ndimage.label(mask)\n",
    "\n",
    "# Find the largest connected component\n",
    "sizes = ndimage.sum(mask, label_im, range(nb_labels + 1))\n",
    "mask_size = sizes < 1000\n",
    "remove_pixel = mask_size[label_im]\n",
    "label_im[remove_pixel] = 0\n",
    "labels = np.unique(label_im)\n",
    "label_im = np.searchsorted(labels, label_im)\n",
    "\n",
    "# Now that we have only one connected component, extract it's bounding box\n",
    "slice_x, slice_y = ndimage.find_objects(label_im==4)[0]\n",
    "roi = im[slice_x, slice_y]\n",
    "\n",
    "plt.figure(figsize=(4, 2))\n",
    "plt.axes([0, 0, 1, 1])\n",
    "plt.imshow(roi)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(900, 900)"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "remove_pixel.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "edges1 = feature.canny(im)\n",
    "edges2 = feature.canny(im, sigma=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<input>:11: DeprecationWarning: invalid escape sequence \\s\n",
      "<input>:11: DeprecationWarning: invalid escape sequence \\s\n",
      "<input>:11: DeprecationWarning: invalid escape sequence \\s\n",
      "<input>:11: DeprecationWarning: invalid escape sequence \\s\n",
      "<input>:11: DeprecationWarning: invalid escape sequence \\s\n",
      "<input>:11: DeprecationWarning: invalid escape sequence \\s\n",
      "<input>:11: DeprecationWarning: invalid escape sequence \\s\n",
      "<ipython-input-22-d4c7fe1e591b>:11: DeprecationWarning: invalid escape sequence \\s\n",
      "  ax2.set_title('Canny filter, $\\sigma=1$', fontsize=20)\n"
     ]
    }
   ],
   "source": [
    "# display results\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(8, 3),\n",
    "                                    sharex=True, sharey=True)\n",
    "\n",
    "ax1.imshow(im, cmap=plt.cm.gray)\n",
    "ax1.axis('off')\n",
    "ax1.set_title('noisy image', fontsize=20)\n",
    "\n",
    "ax2.imshow(edges1, cmap=plt.cm.gray)\n",
    "ax2.axis('off')\n",
    "ax2.set_title('Canny filter, $\\sigma=1$', fontsize=20)\n",
    "\n",
    "ax3.imshow(edges2, cmap=plt.cm.gray)\n",
    "ax3.axis('off')\n",
    "ax3.set_title('Canny filter, $\\sigma=3$', fontsize=20)\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.morphology import convex_hull_image\n",
    "from skimage import data, img_as_float\n",
    "from skimage.util import invert\n",
    "\n",
    "# The original image is inverted as the object must be white.\n",
    "#image = invert(frames[0])\n",
    "image = frames[0]\n",
    "\n",
    "chull = convex_hull_image(image)\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(8, 4))\n",
    "ax = axes.ravel()\n",
    "\n",
    "ax[0].set_title('Original picture')\n",
    "ax[0].imshow(image, cmap=plt.cm.gray, interpolation='nearest')\n",
    "ax[0].set_axis_off()\n",
    "\n",
    "ax[1].set_title('Transformed picture')\n",
    "ax[1].imshow(chull, cmap=plt.cm.gray, interpolation='nearest')\n",
    "ax[1].set_axis_off()\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from skimage.data import camera\n",
    "from skimage.filters import roberts, sobel, scharr, prewitt\n",
    "\n",
    "\n",
    "image = frames[0]\n",
    "edge_roberts = roberts(image)\n",
    "edge_sobel = sobel(image)\n",
    "\n",
    "fig, ax = plt.subplots(ncols=2, sharex=True, sharey=True,\n",
    "                       figsize=(8, 4))\n",
    "\n",
    "ax[0].imshow(edge_roberts, cmap=plt.cm.gray)\n",
    "ax[0].set_title('Roberts Edge Detection')\n",
    "\n",
    "ax[1].imshow(edge_sobel, cmap=plt.cm.gray)\n",
    "ax[1].set_title('Sobel Edge Detection')\n",
    "\n",
    "for a in ax:\n",
    "    a.axis('off')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from skimage.color import rgb2gray\n",
    "from skimage import data\n",
    "from skimage.filters import gaussian\n",
    "from skimage.segmentation import active_contour\n",
    "\n",
    "# Test scipy version, since active contour is only possible\n",
    "# with recent scipy version\n",
    "import scipy\n",
    "split_version = scipy.__version__.split('.')\n",
    "if not(split_version[-1].isdigit()): # Remove dev string if present\n",
    "        split_version.pop()\n",
    "scipy_version = list(map(int, split_version))\n",
    "new_scipy = scipy_version[0] > 0 or \\\n",
    "            (scipy_version[0] == 0 and scipy_version[1] >= 14)\n",
    "\n",
    "img = frames[0]#data.astronaut()\n",
    "img = rgb2gray(img)\n",
    "\n",
    "x = [100, 800, 800, 100, 100]\n",
    "y = [100, 100, 800, 800, 100]\n",
    "init = np.array([x, y]).T\n",
    "\n",
    "if not new_scipy:\n",
    "    print('You are using an old version of scipy. '\n",
    "          'Active contours is implemented for scipy versions '\n",
    "          '0.14.0 and above.')\n",
    "\n",
    "if new_scipy:\n",
    "    snake = active_contour(gaussian(img, 1),\n",
    "                           init, alpha=0.015, beta=10, gamma=0.001)\n",
    "\n",
    "    fig = plt.figure(figsize=(7, 7))\n",
    "    ax = fig.add_subplot(111)\n",
    "    plt.gray()\n",
    "    ax.imshow(img)\n",
    "    ax.plot(init[:, 0], init[:, 1], '--r', lw=3)\n",
    "    ax.plot(snake[:, 0], snake[:, 1], '-b', lw=3)\n",
    "    ax.set_xticks([]), ax.set_yticks([])\n",
    "    ax.axis([0, img.shape[1], img.shape[0], 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0xd236828>"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from skimage import data, color\n",
    "from skimage.transform import hough_circle, hough_circle_peaks\n",
    "from skimage.feature import canny\n",
    "from skimage.draw import circle_perimeter\n",
    "from skimage.util import img_as_ubyte\n",
    "\n",
    "\n",
    "# Load picture and detect edges\n",
    "im = frames[0]\n",
    "image = im#img_as_ubyte(im/im.max())\n",
    "edges = canny(image, sigma=20)#, low_threshold=10, high_threshold=50)\n",
    "\n",
    "\n",
    "# Detect two radii\n",
    "hough_radii = np.arange(5, 100, 5)\n",
    "hough_res = hough_circle(edges, hough_radii)\n",
    "\n",
    "# Select the most prominent 5 circles\n",
    "accums, cx, cy, radii = hough_circle_peaks(hough_res, hough_radii,\n",
    "                                           total_num_peaks=100)\n",
    "\n",
    "# Draw them\n",
    "fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(10, 4))\n",
    "image = color.gray2rgb(image)\n",
    "for center_y, center_x, radius in zip(cy, cx, radii):\n",
    "    circy, circx = circle_perimeter(center_y, center_x, radius)\n",
    "    image[circy, circx] = (220, 20, 20)\n",
    "\n",
    "ax.imshow(image, cmap=plt.cm.gray)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-42-d80450b2c2f6>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     18\u001b[0m \u001b[1;31m# The threshold eliminates low accumulators\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     19\u001b[0m result = hough_ellipse(edges, accuracy=20, threshold=250,\n\u001b[1;32m---> 20\u001b[1;33m                        min_size=100, max_size=120)\n\u001b[0m\u001b[0;32m     21\u001b[0m \u001b[0mresult\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msort\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0morder\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'accumulator'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     22\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Anaconda3\\envs\\wradlib0_11\\lib\\site-packages\\skimage\\transform\\hough_transform.py\u001b[0m in \u001b[0;36mhough_ellipse\u001b[1;34m(img, threshold, accuracy, min_size, max_size)\u001b[0m\n\u001b[0;32m    161\u001b[0m     \"\"\"\n\u001b[0;32m    162\u001b[0m     return _hough_ellipse(img, threshold=threshold, accuracy=accuracy,\n\u001b[1;32m--> 163\u001b[1;33m                           min_size=min_size, max_size=max_size)\n\u001b[0m\u001b[0;32m    164\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    165\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mskimage/transform/_hough_transform.pyx\u001b[0m in \u001b[0;36mskimage.transform._hough_transform._hough_ellipse (skimage\\transform\\_hough_transform.c:4430)\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;32mC:\\Anaconda3\\envs\\wradlib0_11\\lib\\site-packages\\numpy\\core\\fromnumeric.py\u001b[0m in \u001b[0;36mamax\u001b[1;34m(a, axis, out, keepdims)\u001b[0m\n\u001b[0;32m   2305\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2306\u001b[0m     \"\"\"\n\u001b[1;32m-> 2307\u001b[1;33m     \u001b[0mkwargs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m{\u001b[0m\u001b[1;33m}\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   2308\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0mkeepdims\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_NoValue\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2309\u001b[0m         \u001b[0mkwargs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'keepdims'\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from skimage import data, color, img_as_ubyte\n",
    "from skimage.feature import canny\n",
    "from skimage.transform import hough_ellipse\n",
    "from skimage.draw import ellipse_perimeter\n",
    "\n",
    "# Load picture, convert to grayscale and detect edges\n",
    "#image_rgb = data.coffee()[0:220, 160:420]\n",
    "#image_gray = color.rgb2gray(image_rgb)\n",
    "image_gray = frames[0]\n",
    "edges = canny(image_gray, sigma=2.0,\n",
    "              low_threshold=0.55, high_threshold=0.8)\n",
    "\n",
    "# Perform a Hough Transform\n",
    "# The accuracy corresponds to the bin size of a major axis.\n",
    "# The value is chosen in order to get a single high accumulator.\n",
    "# The threshold eliminates low accumulators\n",
    "result = hough_ellipse(edges, accuracy=20, threshold=250,\n",
    "                       min_size=100, max_size=120)\n",
    "result.sort(order='accumulator')\n",
    "\n",
    "# Estimated parameters for the ellipse\n",
    "best = list(result[-1])\n",
    "yc, xc, a, b = [int(round(x)) for x in best[1:5]]\n",
    "orientation = best[5]\n",
    "\n",
    "# Draw the ellipse on the original image\n",
    "cy, cx = ellipse_perimeter(yc, xc, a, b, orientation)\n",
    "image_rgb[cy, cx] = (0, 0, 255)\n",
    "# Draw the edge (white) and the resulting ellipse (red)\n",
    "edges = color.gray2rgb(img_as_ubyte(edges))\n",
    "edges[cy, cx] = (250, 0, 0)\n",
    "\n",
    "fig2, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(8, 4), sharex=True,\n",
    "                                sharey=True,\n",
    "                                subplot_kw={'adjustable':'box-forced'})\n",
    "\n",
    "ax1.set_title('Original picture')\n",
    "ax1.imshow(image_rgb)\n",
    "\n",
    "ax2.set_title('Edge (white) and result (red)')\n",
    "ax2.imshow(edges)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'25 DAISY descriptors extracted:')"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from skimage.feature import daisy\n",
    "from skimage import data\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "img = frames[0]\n",
    "descs, descs_img = daisy(img, step=180, radius=58, rings=2, histograms=6,\n",
    "                         orientations=8, visualize=True)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.axis('off')\n",
    "ax.imshow(descs_img)\n",
    "descs_num = descs.shape[0] * descs.shape[1]\n",
    "ax.set_title('%i DAISY descriptors extracted:' % descs_num)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use OpenCV's Optical Flow to detect and track features\n",
    "\n",
    "This example uses the Lucas-Kanade Optical Flow implementation in OpenCV (see [here](http://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_video/py_lucas_kanade/py_lucas_kanade.html)). We take the first frame, detect some Shi-Tomasi corner points in it, then we iteratively track those points over the subsequent images.\n",
    "\n",
    "The parameter dictionaries are certainly something to experiment with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# FEATURE DETECTION: Parameters for ShiTomasi corner detection\n",
    "feature_params = dict( maxCorners = 200,\n",
    "                       qualityLevel = 0.2,\n",
    "                       minDistance = 7,\n",
    "                       blockSize = 21 )\n",
    "\n",
    "# FEATURE TRACKING: Parameters for Lucas Kanade (lk) Optical Flow technique\n",
    "lk_params = dict( winSize  = (20,20),\n",
    "                  maxLevel = 2,\n",
    "                  criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 10, 0))\n",
    "\n",
    "# Over which time steps (of the data we've read in) do you want to track\n",
    "trackstart = 0\n",
    "trackend = 36"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Our approach requires 8 bit integers - so we need to normalize our radar data accordingly\n",
    "#   (there might be a more elegant solution...)\n",
    "minval = 0\n",
    "maxval = 59 # dBZ in this case\n",
    "iframes = frames.copy()\n",
    "iframes[iframes<minval] = minval\n",
    "iframes[iframes>maxval] = maxval\n",
    "iframes = ((iframes / maxval)*255).astype(np.uint8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find good features to track...\n",
    "old = cv2.goodFeaturesToTrack(iframes[trackstart], mask = None, **feature_params)\n",
    "print(\"Found %d good features to track.\" % len(old) )\n",
    "\n",
    "# Set containers to collect results (time steps in rows, detected corners in columns)\n",
    "#   Tracking status\n",
    "sts = np.zeros((trackend,len(old)), dtype=np.bool)\n",
    "#   corner x coords\n",
    "x = np.zeros((trackend,len(old))) * np.nan\n",
    "#   corner y coords\n",
    "y = np.zeros((trackend,len(old))) * np.nan\n",
    "#   tracking error\n",
    "errs = np.zeros((trackend,len(old))) * np.nan\n",
    "#   Assign persistent corner IDs\n",
    "ids = np.arange(len(old))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Track good features\n",
    "for i in range(trackstart, trackend):\n",
    "    # track current corners in next image\n",
    "    new, st, err = cv2.calcOpticalFlowPyrLK(prevImg=iframes[i], nextImg=iframes[i+1], prevPts=old, nextPts=None, **lk_params)\n",
    "    success = st.ravel()==1\n",
    "    ids = ids[success]\n",
    "    sts[i, ids] = True\n",
    "    x[i, ids] = old[success,0,0]\n",
    "    y[i, ids] = old[success,0,1]\n",
    "    errs[i, ids] = err.ravel()[success]\n",
    "    # new corners will be old in the next loop\n",
    "    old = new[success]\n",
    "\n",
    "# Incremental euclidic distance from starting point\n",
    "trackdist = np.diff( np.sqrt( (x-x[0].reshape((1,-1)))**2 + (y-y[0].reshape((1,-1)))**2 ), axis=0 )\n",
    "trackdist = np.vstack( (np.zeros((1,trackdist.shape[1])), trackdist))\n",
    "\n",
    "# Plot feature persistence\n",
    "fig = plt.figure(figsize=(12,6))\n",
    "ax = fig.add_subplot(211)\n",
    "cb = plt.imshow(errs, interpolation=\"none\", cmap=\"summer\", vmax = 15)\n",
    "plt.xlabel(\"Feature ID\")\n",
    "plt.ylabel(\"Tracking time step\")\n",
    "plt.colorbar(cb, shrink=0.5)\n",
    "plt.title(\"Tracking error\")\n",
    "\n",
    "# Plot consistence of movement\n",
    "ax = fig.add_subplot(212)\n",
    "cb = plt.imshow(trackdist, interpolation=\"none\", cmap=\"bwr\", vmin=-5, vmax=5)\n",
    "plt.xlabel(\"Feature ID\")\n",
    "plt.ylabel(\"Tracking time step\")\n",
    "plt.colorbar(cb, shrink=0.75)\n",
    "plt.title(\"Incremental euclidian distance from starting point\")\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find good tracks (but what is a \"good\" track...?)\n",
    "#   Certainly a lot of subjective criteria to play with...\n",
    "goodtrack = np.zeros(x.shape[1], dtype=np.bool)\n",
    "for i in range(len(goodtrack)):\n",
    "    # persistence of the track\n",
    "    if len(np.where(sts[:,i])[0]) < 2:\n",
    "        continue\n",
    "    # consistency of movement\n",
    "    if len(np.where(trackdist[:,i]<0)[0]) > 0:\n",
    "        continue\n",
    "    # tracking error\n",
    "    if len(np.where(errs[:,i]>15)[0]) > 5:\n",
    "        continue\n",
    "    goodtrack[i] = True\n",
    "print(\"Found %d good tracks and %d bad tracks.\" % \\\n",
    "      (len(np.where(goodtrack)[0]), len(goodtrack)-len(np.where(goodtrack)[0])) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize tracks: green=good track, red=bad track\n",
    "goodcolor = \"limegreen\"\n",
    "badcolor = \"red\"\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "ax = fig.add_subplot(111, aspect=\"equal\")\n",
    "# average reflectivity over entire tracking period as background image\n",
    "ax.imshow(np.mean(frames[trackstart:trackend], axis=0), origin=\"lower\", cmap=\"gray\", interpolation=\"none\")\n",
    "plt.xlabel(\"Easting (# pixels)\")\n",
    "plt.ylabel(\"Northing (# pixels)\")\n",
    "plt.title(\"[Zoom in to inspect track properties (not in inline mode!)]\")\n",
    "plt.grid(color=\"white\")\n",
    "plt.xlim(0,nx/downsizeby)\n",
    "plt.ylim(0,nx/downsizeby)\n",
    "bad_line = plt.Line2D([], [], color=badcolor, label='Bad track')\n",
    "good_line = plt.Line2D([], [], color=goodcolor, label='Good track')\n",
    "plt.legend(handles=[bad_line, good_line], loc=\"upper left\")\n",
    "for i, isgood in enumerate(goodtrack):\n",
    "    ix = sts[:,i]\n",
    "    color = badcolor\n",
    "    if isgood:\n",
    "        color = goodcolor\n",
    "    ax.plot(x[ix,i], y[ix,i],marker=\"None\", color=color, markersize=14, linestyle=\"-\")\n",
    "    ax.arrow(x[ix,i][-2], y[ix,i][-2],\n",
    "             np.diff(x[ix,i][-2:])[0], np.diff(y[ix,i][-2:])[0], \n",
    "             head_width=2, head_length=2, fc=color, ec=color)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Animate features\n",
    "\n",
    "# Prepare canvas\n",
    "fig = plt.figure(figsize=(8,8))\n",
    "ax1 = plt.subplot(111,aspect=\"equal\")\n",
    "im1 = ax1.imshow(iframes[trackstart], origin=\"lower\", cmap=\"gray\", interpolation=\"none\")\n",
    "plt.xlabel(\"Easting (# pixels)\")\n",
    "plt.ylabel(\"Northing (# pixels)\")\n",
    "plt.title(\"[Zoom in to inspect track properties (not in inline mode!)]\")\n",
    "plt.grid(color=\"white\")\n",
    "plt.xlim(0,nx/downsizeby)\n",
    "plt.ylim(0,nx/downsizeby)\n",
    "ax1.plot(x[0,goodtrack], y[0,goodtrack], linestyle=\"None\", marker=\"o\", mfc=\"None\", mec=\"limegreen\")\n",
    "ax1.plot(x[0,~goodtrack], y[0,~goodtrack], linestyle=\"None\", marker=\"o\", mfc=\"None\", mec=\"red\")\n",
    "ax1.grid(color=\"white\")\n",
    "tstamp1 = ax1.text(25, 850, dtimes[trackstart].isoformat(), color=\"white\", fontsize=14)\n",
    "\n",
    "def animate(j):\n",
    "    im1.set_array(iframes[trackstart+j])\n",
    "    for line in plt.gca().get_lines():\n",
    "        if not line.get_linestyle()==\"None\":\n",
    "            line.remove()\n",
    "    for i, isgood in enumerate(goodtrack):\n",
    "        ix = np.where(sts[:j,i])[0]\n",
    "        color = \"red\"\n",
    "        if isgood:\n",
    "            color = \"limegreen\"\n",
    "        ax1.plot(x[ix,i], y[ix,i], marker=\"None\", color=color, markersize=14, linestyle=\"-\")\n",
    "    tstamp1.set_text(dtimes[trackstart+j].isoformat())\n",
    "    return im1\n",
    "\n",
    "# ATTENTION: THIS IS SLOW - Rendering each frame of the animation might take more time than the interval between the frames\n",
    "#    This can cause the temporal sequence to be confused in the matplotlib interactive mode.\n",
    "#    The animation thus looks better if saved as movie, or you have to increase the interval argument\n",
    "#    Animation not shown in notebook if you use %pylab inline\n",
    "ani = animation.FuncAnimation(fig, animate, frames=np.arange(trackstart, trackend-1), interval=400, blit=False)\n",
    "#ani.save(\"features.avi\", dpi=500, bitrate=2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Update tracked corners for each time step of the considered tracking period\n",
    "\n",
    "Until now, we only tracked those corners which we detected in the initial time step. We now want to add new tracks with each addtional time step, and follow these as well. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_crns = [cv2.goodFeaturesToTrack(iframes[i], mask = None, **feature_params) for i in range(trackstart, trackend)]\n",
    "print(\"List of # corners in each time step:\\n\", [len(crn) for crn in init_crns ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# this function wraps up everything which we already did above for a single set of corners\n",
    "def tracker(old, frameset, lk_params):\n",
    "    # Set containers to collect results (time steps in rows, corners in columns)\n",
    "    #   Tracking status\n",
    "    sts = np.zeros((trackend,len(old)), dtype=np.bool)\n",
    "    #   corner x coords\n",
    "    x = np.zeros((trackend,len(old))) * np.nan\n",
    "    #   corner y coords\n",
    "    y = np.zeros((trackend,len(old))) * np.nan\n",
    "    #   tracking error\n",
    "    errs = np.zeros((trackend,len(old))) * np.nan\n",
    "    #   Assign persistent corner IDs\n",
    "    ids = np.arange(len(old))\n",
    "    # Track good features\n",
    "    for i in range(len(frameset)-1):\n",
    "        # track current corners in next image\n",
    "        new, st, err = cv2.calcOpticalFlowPyrLK(prevImg=frameset[i], nextImg=frameset[i+1],\n",
    "                                                prevPts=old, nextPts=None, **lk_params)\n",
    "        success = st.ravel()==1\n",
    "        ids = ids[success]\n",
    "        sts[i, ids] = True\n",
    "        x[i, ids] = new[success,0,0]\n",
    "        y[i, ids] = new[success,0,1]\n",
    "        errs[i, ids] = err.ravel()[success]\n",
    "        # new corners will be old in the next loop\n",
    "        old = new[success]\n",
    "\n",
    "    # Incremental euclidic distance from starting point\n",
    "    trackdist = np.diff( np.sqrt( (x-x[0].reshape((1,-1)))**2 + (y-y[0].reshape((1,-1)))**2 ), axis=0 )\n",
    "    trackdist = np.vstack( (np.zeros((1,trackdist.shape[1])), trackdist))\n",
    "\n",
    "    # Find good tracks (but what is a \"good\" track...?)\n",
    "    goodtrack = np.zeros(x.shape[1], dtype=np.bool)\n",
    "    for i in range(len(goodtrack)):\n",
    "        # persistence of the track\n",
    "        if len(np.where(sts[:,i])[0]) < 2:\n",
    "            continue\n",
    "        # consistency of movement\n",
    "        if len(np.where(trackdist[:,i]<0)[0]) > 0:\n",
    "            continue\n",
    "        # tracking error\n",
    "        if len(np.where(errs[:,i]>15)[0]) > 5:\n",
    "            continue\n",
    "        goodtrack[i] = True\n",
    "    \n",
    "    return sts, x, y, errs, goodtrack    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sts_ls, x_ls, y_ls, errs_ls, goodtrack_ls = [], [], [], [], []\n",
    "for i, crns in enumerate(init_crns):\n",
    "    sts, x, y, errs, goodtrack = tracker(crns, iframes[i:], lk_params)\n",
    "    sts_ls.append(sts)\n",
    "    x_ls.append(x)\n",
    "    y_ls.append(y)\n",
    "    errs_ls.append(errs)\n",
    "    goodtrack_ls.append(goodtrack)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize tracks:\n",
    "fig = plt.figure(figsize=(12,12))\n",
    "ax = fig.add_subplot(111, aspect=\"equal\")\n",
    "# average reflectivity as background image\n",
    "ax.imshow(np.mean(frames[trackstart:trackend], axis=0), origin=\"lower\", cmap=\"gray\", interpolation=\"none\")\n",
    "plt.xlabel(\"Easting (# pixels)\")\n",
    "plt.ylabel(\"Northing (# pixels)\")\n",
    "plt.title(\"[Zoom in to inspect track properties (not in inline mode!)]\")\n",
    "plt.grid(color=\"white\")\n",
    "plt.xlim(0,nx/downsizeby)\n",
    "plt.ylim(0,nx/downsizeby)\n",
    "colors = [ plt.cm.spring(i) for i in np.linspace(0,254, len(goodtrack_ls)).astype(\"i4\") ]\n",
    "for j, goodtrack in enumerate(goodtrack_ls[:-2]):\n",
    "    sts, x, y = sts_ls[j], x_ls[j], y_ls[j]\n",
    "    for i, isgood in enumerate(goodtrack):\n",
    "        ix = sts[:,i]\n",
    "        # HERE WE DO NOT PLOT THE BAD TRACKS\n",
    "        color = \"none\"\n",
    "        if isgood:\n",
    "            color = colors[j]\n",
    "        ax.plot(x[ix,i], y[ix,i],marker=\"None\", color=color, linestyle=\"-\", alpha=0.4)\n",
    "        ax.arrow(x[ix,i][-2], y[ix,i][-2],\n",
    "                 np.diff(x[ix,i][-2:])[0], np.diff(y[ix,i][-2:])[0], \n",
    "                 head_width=2, head_length=2, fc=color, ec=color, alpha=0.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ATTENTION: THIS ANIMATION TAKES A LONG WHILE (SEVERAL MINUTES) AND MIGHT STILL BE BUGGY\n",
    "\n",
    "# Prepare canvas\n",
    "fig = plt.figure(figsize=(8,8))\n",
    "ax1 = plt.subplot(111,aspect=\"equal\")\n",
    "im1 = ax1.imshow(iframes[trackstart], origin=\"lower\", cmap=\"gray\", interpolation=\"none\")\n",
    "plt.xlabel(\"Easting (# pixels)\")\n",
    "plt.ylabel(\"Northing (# pixels)\")\n",
    "plt.title(\"[Zoom in to inspect track properties (not in inline mode!)]\")\n",
    "plt.grid(color=\"white\")\n",
    "plt.xlim(0,nx/downsizeby)\n",
    "plt.ylim(0,nx/downsizeby)\n",
    "#ax1.plot(x[0,goodtrack], y[0,goodtrack], linestyle=\"None\", marker=\"o\", mfc=\"None\", mec=colors[0])\n",
    "ax1.grid(color=\"white\")\n",
    "tstamp1 = ax1.text(25, 850, dtimes[trackstart].isoformat(), color=\"white\", fontsize=14)\n",
    "\n",
    "def animate(j):\n",
    "    im1.set_array(iframes[trackstart+j])\n",
    "    for line in plt.gca().get_lines():\n",
    "        line.remove()\n",
    "        #if not line.get_linestyle()==\"None\":\n",
    "        #    line.remove()\n",
    "    for k, goodtrack in enumerate(goodtrack_ls[:j]):\n",
    "        sts, x, y = sts_ls[k], x_ls[k], y_ls[k]\n",
    "        for i, isgood in enumerate(goodtrack):\n",
    "            ix = np.where(sts[:j,i])[0]\n",
    "            # HERE WE DO NOT PLOT THE BAD TRACKS\n",
    "            color = \"none\"\n",
    "            if isgood:\n",
    "                color = colors[k]\n",
    "            #ax1.plot(x[0,goodtrack], y[0,goodtrack], linestyle=\"None\", marker=\"o\", mfc=\"None\", mec=color, alpha=0.4)\n",
    "            ax1.plot(x[ix,i], y[ix,i],marker=\"None\", color=color, linestyle=\"-\", alpha=0.4)\n",
    "# ATTENTION: THIS IS SLOW - Rendering each frame of the animation might take more time than the interval between the frames\n",
    "#    This can cause the temporal sequence to be confused.\n",
    "#    The animation thus looks better if saved as movie, or you have to increase the interval argument\n",
    "ani = animation.FuncAnimation(fig, animate, frames=np.arange(trackstart, trackend-1), interval=200, blit=False)\n",
    "#ani.save(\"features2.avi\", dpi=500, bitrate=2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Experiment with SIFT/SURF feature detection and description\n",
    "\n",
    "See [SIFT and SURF](http://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_feature2d/py_table_of_contents_feature2d/py_table_of_contents_feature2d.html) for feature detection. Right now, this does not seem to add value as compared to the Optical Flow approach above. Features seem to be much less persistent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SURF\n",
    "surf = cv2.xfeatures2d.SURF_create(3000)\n",
    "\n",
    "kplist = []\n",
    "deslist= []\n",
    "\n",
    "for i in range(trackstart, trackend):\n",
    "    kp, des = surf.detectAndCompute(iframes[i],None)\n",
    "    kplist.append(kp)\n",
    "    deslist.append(des)\n",
    "    print(\"Found %d keypoints in step %d.\" % (len(kp), i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10,10))\n",
    "ax = fig.add_subplot(111, aspect=\"equal\")\n",
    "# average reflectivity as background image\n",
    "ax.imshow(frames[0], origin=\"lower\", cmap=\"gray\", interpolation=\"none\")\n",
    "plt.xlabel(\"Easting (# pixels)\")\n",
    "plt.ylabel(\"Northing (# pixels)\")\n",
    "plt.title(\"[Zoom in to inspect feature properties (not in inline mode)]\")\n",
    "tstamp1 = ax1.text(25, 850, dtimes[0].isoformat(), color=\"white\", fontsize=14)\n",
    "plt.grid(color=\"white\")\n",
    "plt.xlim(0,nx/downsizeby)\n",
    "plt.ylim(0,nx/downsizeby)\n",
    "patches = []\n",
    "for kp_ in kplist[0]:\n",
    "    if kp_.size > 5:\n",
    "        circle = mpatches.Circle(kp_.pt, kp_.size, fill=False, edgecolor=\"red\")\n",
    "    #ax.add_patch(circle)\n",
    "    patches.append(circle)\n",
    "collection = PatchCollection(patches, facecolor=\"none\", edgecolor=\"red\")\n",
    "ax.add_collection(collection)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Make list of patch collections for all timesteps\n",
    "def collect(kp):\n",
    "    patches = []\n",
    "    for kp_ in kp:\n",
    "        if (kp_.size > 10) and (kp_.size < 50):\n",
    "            circle = mpatches.Circle(kp_.pt, kp_.size, fill=False, edgecolor=\"red\")\n",
    "            patches.append(circle)\n",
    "    return(PatchCollection(patches, facecolor=\"none\", edgecolor=\"red\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Animate features\n",
    "_plot_style = dict(markersize=12, markeredgewidth=2,\n",
    "                       markerfacecolor='none', markeredgecolor='r',\n",
    "                       marker='o', linestyle='none')\n",
    "_pcm_style = dict(cmap=plt.cm.spectral, vmin=0., vmax=30.)\n",
    "\n",
    "# Prepare canvas\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "ax1 = plt.subplot(111,aspect=\"equal\")\n",
    "im1 = ax1.imshow(iframes[0], origin=\"lower\", cmap=\"gray\", interpolation=\"none\")\n",
    "ax1.add_collection(collect(kplist[0]))\n",
    "ax1.grid(color=\"white\")\n",
    "tstamp1 = ax1.text(25, 850, dtimes[0].isoformat(), color=\"white\", fontsize=14)\n",
    "\n",
    "\n",
    "def animate(i):\n",
    "    im1.set_array(iframes[trackstart+i])\n",
    "    ax1.collections = []\n",
    "    ax1.add_collection(collect(kplist[trackstart+i]))\n",
    "    tstamp1.set_text(dtimes[trackstart+i].isoformat())\n",
    "    return im1\n",
    "\n",
    "ani = animation.FuncAnimation(fig, animate, frames=np.arange(trackstart, trackend-1), interval=200, blit=False)\n",
    "#ani.save(\"features_surf.avi\", dpi=400, bitrate=2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Match features (brute force)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According [Bowler et al. (2004)](http://www.sciencedirect.com/science/article/pii/S0022169403004591), maximum advection velocity of rainfall objects is about 130 km/h which is roughly 10 km (pixels) in 5 minutes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "maxveloc = 10.\n",
    "# Detect initial feature set \n",
    "detector = cv2.xfeatures2d.SURF_create(3000)\n",
    "kp1, des1 = detector.detectAndCompute(iframes[trackstart],None)\n",
    "\n",
    "# create BFMatcher object\n",
    "bf = cv2.BFMatcher()\n",
    "\n",
    "kp1_ls = []\n",
    "kp2_ls = []\n",
    "\n",
    "for i in range(trackstart+1, trackend):\n",
    "    kp2, des2 = detector.detectAndCompute(iframes[i],None)\n",
    "    matches = bf.knnMatch(des1, des2, k=1)\n",
    "    # Select matches to keep\n",
    "    kp1_, des1_, kp2_, des2_  = [], [], [], []\n",
    "    for match in matches:\n",
    "        match=match[0]\n",
    "        xy = np.vstack( (kp1[match.queryIdx].pt, kp2[match.trainIdx].pt) )\n",
    "        eucdist = np.sqrt( (xy[0,0] - xy[1,0])**2 + (xy[0,1] - xy[1,1])**2 )\n",
    "        if eucdist < maxveloc:\n",
    "            kp1_.append( kp1[match.queryIdx] )\n",
    "            des1_.append( np.array( des1[match.queryIdx] ) )\n",
    "            kp2_.append( kp2[match.trainIdx] )\n",
    "            des2_.append( np.array( des2[match.trainIdx] ) )\n",
    "    kp1_ls.append(kp1_)\n",
    "    kp2_ls.append(kp2_)\n",
    "    # Update initial feature set\n",
    "    kp1, des1 = kp2_, np.array( des2_ )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.4"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "119px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}