{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "06786cfc-624c-d137-3496-61cb37023558"
      },
      "source": [
        "# Intro\n",
        "Step-by-step guide for extracting features from shapes by turning them into time-series. The functions are optimized for the Swedish Leaf Dataset as it is published on Kaggle.\n",
        "\n",
        "Data Scientists spend vast majority of their time by data preparation, not model optimization. So let's dive into that a bit.\n",
        "\n",
        "## [Here is the second part][1], it is more advanced and [follows PEP8][2]\n",
        "\n",
        "![Feature Extraction from leaf contours][3]\n",
        "\n",
        "\n",
        "\n",
        "  [1]: https://www.kaggle.com/lorinc/leaf-classification/feature-extraction-v4/\n",
        "  [2]: https://www.python.org/dev/peps/pep-0008/\n",
        "  [3]: https://raw.githubusercontent.com/lorinc/kaggle-notebooks/master/extracted_leaf_shape.png"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "6e96dca3-b41b-5410-edf3-544656e86067"
      },
      "outputs": [],
      "source": [
        "# imports\n",
        "import numpy as np                     # numeric python lib\n",
        "\n",
        "import matplotlib.image as mpimg       # reading images to numpy arrays\n",
        "import matplotlib.pyplot as plt        # to plot any graph\n",
        "import matplotlib.patches as mpatches  # to draw a circle at the mean contour\n",
        "\n",
        "from skimage import measure            # to find shape contour\n",
        "import scipy.ndimage as ndi            # to determine shape centrality\n",
        "\n",
        "\n",
        "# matplotlib setup\n",
        "%matplotlib inline\n",
        "from pylab import rcParams\n",
        "rcParams['figure.figsize'] = (6, 6)      # setting default size of plots"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "4e58ed89-aad1-955a-92a5-5ddcf91584a4"
      },
      "outputs": [],
      "source": [
        "# reading an image file using matplotlib into a numpy array\n",
        "# good ones: 11, 19, 23, 27, 48, 53, 78, 218\n",
        "img = mpimg.imread('../input/images/53.jpg')\n",
        "\n",
        "# using image processing module of scipy to find the center of the leaf\n",
        "cy, cx = ndi.center_of_mass(img)\n",
        "\n",
        "plt.imshow(img, cmap='Set3')  # show me the leaf\n",
        "plt.scatter(cx, cy)           # show me its center\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "e80c287f-d874-571b-c299-b474668c0ffb"
      },
      "source": [
        "Later we might want to switch to another measure of centrality, based on how efficient this center is, when we generate a time-series from the shape, using the distance between the edge and the center.\n",
        "\n",
        "One way to do that is just measure the (Euclidean) distance between the center and the edge... but there is a better way - we project the Cartesian coordinates into Polar coordinates.\n",
        "\n",
        "But before that, we need to [find the edges][1] of the leaf.\n",
        "\n",
        "  [1]: http://scikit-image.org/docs/dev/auto_examples/"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "294e98a4-1766-508b-4ba9-89e26c1ae8d5"
      },
      "outputs": [],
      "source": [
        "# scikit-learn imaging contour finding, returns a list of found edges\n",
        "contours = measure.find_contours(img, .8)\n",
        "\n",
        "# from which we choose the longest one\n",
        "contour = max(contours, key=len)\n",
        "\n",
        "# let us see the contour that we hopefully found\n",
        "plt.plot(contour[::,1], contour[::,0], linewidth=0.5)  # (I will explain this [::,x] later)\n",
        "plt.imshow(img, cmap='Set3')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "943605c1-46d8-fedf-da8a-7e55525b8559"
      },
      "source": [
        "Now we can project this contour - that is, thousands of pairs of (x,y) coordinates - into the Polar coordinate system.\n",
        "![Cartesian to Polar][1]\n",
        "[1]: http://jwilson.coe.uga.edu/emat6680fa11/lee/asnmt11hylee/fig2.jpg"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "863a55db-818b-b1e8-5bee-41293a5b796e"
      },
      "outputs": [],
      "source": [
        "# cartesian to polar coordinates, just as the image shows above\n",
        "def cart2pol(x, y):\n",
        "    rho = np.sqrt(x**2 + y**2)\n",
        "    phi = np.arctan2(y, x)\n",
        "    return [rho, phi]\n",
        "\n",
        "# just calling the transformation on all pairs in the set\n",
        "polar_contour = np.array([cart2pol(x, y) for x, y in contour])\n",
        "\n",
        "# and plotting the result\n",
        "plt.plot(polar_contour[::,1], polar_contour[::,0], linewidth=0.5)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "52967ca3-dd17-7256-9008-e6c332539a78"
      },
      "source": [
        "What on earth is this. We expected a time-series, then we got a monster-leaf.\n",
        "\n",
        "Let's try the projection again, but move the leaf to (0,0) first this time.\n",
        "\n",
        "This gives a chance to look into something we already used, but not explained: [numpy.ndarray indexing][1]. \"nd\" stands for \"n-dimensional\" and there is a powerful syntax helps us to select whatever slice we need from our array. The syntax is `array[start:stop:step]`. Also, this syntax allows us to select into subdimensions.\n",
        "\n",
        "\n",
        "  [1]: http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "1b0f8711-3150-6f44-ecee-dc2c36568d5d"
      },
      "outputs": [],
      "source": [
        "# numpy BASIC indexing example, see link above for more\n",
        "\n",
        "x = np.array([[[1,11,111], [2,22,222], [3,33,333]], \n",
        "              [[4,44,444], [5,55,555], [6,66,666]], \n",
        "              [[7,77,777], [8,88,888], [9,99,999]]])\n",
        "\n",
        "# reverse the first dimension\n",
        "# take the 0th element\n",
        "# and take its last element\n",
        "x[::-1, 0, -1]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "69b72c07-ed38-0d37-e2c3-dffd9e54359e"
      },
      "source": [
        "So, using this, let us demean the contour data.\n",
        "\n",
        "Demean (making its mean 0 by shifting points) is needed: the polar coordnate projection failed to yield what we want, because the shape is in the +,+ part of the Cartesian system, not around the center.\n",
        "\n",
        "Why the contour, why not the image? First, image is big, contour is small. Second, image is just hundreds n hundreds of bits on a grid, not (x,y) pairs. It is not the right data format for us."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "ef1efde5-32aa-c259-710d-ff95cf483661"
      },
      "outputs": [],
      "source": [
        "# numpy is smart and assumes the same about us\n",
        "# if we substract a number from an array of numbers,\n",
        "# it assumes that we wanted to substract from all members\n",
        "contour[::,1] -= cx  # demean X\n",
        "contour[::,0] -= cy  # demean Y"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "303e4a6a-c9a5-c498-15e3-3080f39bb2d3"
      },
      "outputs": [],
      "source": [
        "# checking if we succeeded to move the center to (0,0)\n",
        "plt.plot(-contour[::,1], -contour[::,0], linewidth=0.5)\n",
        "plt.grid()\n",
        "plt.scatter(0, 0)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "d1013727-94ae-9cc3-1262-65ed20beaf40"
      },
      "source": [
        "Now we can try to project it again into polar space."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "81cdbf38-207c-5575-83df-1e5a77032931"
      },
      "outputs": [],
      "source": [
        "# just calling the transformation on all pairs in the set\n",
        "polar_contour = np.array([cart2pol(x, y) for x, y in contour])\n",
        "\n",
        "# and plotting the result\n",
        "rcParams['figure.figsize'] = (12, 6)\n",
        "plt.subplot(121)\n",
        "plt.scatter(polar_contour[::,1], polar_contour[::,0], linewidth=0, s=.5, c=polar_contour[::,1])\n",
        "plt.title('in Polar Coordinates')\n",
        "plt.grid()\n",
        "plt.subplot(122)\n",
        "plt.scatter(contour[::,1],             # x axis is radians\n",
        "            contour[::,0],             # y axis is distance from center\n",
        "            linewidth=0, s=2,          # small points, w/o borders\n",
        "            c=range(len(contour)))     # continuous coloring (so that plots match)\n",
        "plt.scatter(0, 0)\n",
        "plt.title('in Cartesian Coordinates')\n",
        "plt.grid()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "1a1c2893-da6a-57fe-316f-c0cbdd01903e"
      },
      "source": [
        "Ok, I consider this beginning part done.\n",
        "\n",
        "Challenges:\n",
        "\n",
        "- This is definitely not a time-series yet, as one x can have multiple y values. (I do not know yet if it is a good thing for the classification and feature extraction or not)\n",
        "- I have to solve that the same leaf if I rotate it, is still equivalent from feature extraction and classification point of view.\n",
        "- And - of course - I'm yet to find and quantify meaningful and generic features.\n",
        "\n",
        "And a funny thought - what if I cut each sample into half along the symmetry axis and treat them as two samples? How shameless would that be? I get double sample size, I generalize better, I simplify and reduce the features I work with... \n",
        "\n",
        "![Evil genius laughter][1]\n",
        "\n",
        "  [1]: http://www.jimbowley.com/wp-content/uploads/2011/06/DrEvilPinky.jpg"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "a7c8f072-4cb1-8f7a-c3af-63b124b172d9"
      },
      "outputs": [],
      "source": [
        "# check a few scikitlearn image feature extractions, if they can help us\n",
        "\n",
        "from skimage.feature import corner_harris, corner_subpix, corner_peaks, CENSURE\n",
        "\n",
        "detector = CENSURE()\n",
        "detector.detect(img)\n",
        "\n",
        "coords = corner_peaks(corner_harris(img), min_distance=5)\n",
        "coords_subpix = corner_subpix(img, coords, window_size=13)\n",
        "\n",
        "plt.subplot(121)\n",
        "plt.title('CENSURE feature detection')\n",
        "plt.imshow(img, cmap='Set3')\n",
        "plt.scatter(detector.keypoints[:, 1], detector.keypoints[:, 0],\n",
        "              2 ** detector.scales, facecolors='none', edgecolors='r')\n",
        "\n",
        "plt.subplot(122)\n",
        "plt.title('Harris Corner Detection')\n",
        "plt.imshow(img, cmap='Set3')  # show me the leaf\n",
        "plt.plot(coords[:, 1], coords[:, 0], '.b', markersize=5)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "f4428744-e69b-dd6c-bffd-dda5d69212b0"
      },
      "source": [
        "I really want to find features that generalize well, but this is just random noise.\n",
        "\n",
        "How about finding local maxima and minima on the time series instead?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "2072c0c8-eead-b399-4703-710d15367af6"
      },
      "outputs": [],
      "source": [
        "from scipy.signal import argrelextrema\n",
        "\n",
        "# for local maxima\n",
        "c_max_index = argrelextrema(polar_contour[::,0], np.greater, order=50)\n",
        "c_min_index = argrelextrema(polar_contour[::,0], np.less, order=50)\n",
        "\n",
        "plt.subplot(121)\n",
        "plt.scatter(polar_contour[::,1], polar_contour[::,0], \n",
        "            linewidth=0, s=2, c='k')\n",
        "plt.scatter(polar_contour[::,1][c_max_index], \n",
        "            polar_contour[::,0][c_max_index], \n",
        "            linewidth=0, s=30, c='b')\n",
        "plt.scatter(polar_contour[::,1][c_min_index], \n",
        "            polar_contour[::,0][c_min_index], \n",
        "            linewidth=0, s=30, c='r')\n",
        "\n",
        "plt.subplot(122)\n",
        "plt.scatter(contour[::,1], contour[::,0], \n",
        "            linewidth=0, s=2, c='k')\n",
        "plt.scatter(contour[::,1][c_max_index], \n",
        "            contour[::,0][c_max_index], \n",
        "            linewidth=0, s=30, c='b')\n",
        "plt.scatter(contour[::,1][c_min_index], \n",
        "            contour[::,0][c_min_index], \n",
        "            linewidth=0, s=30, c='r')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "1a41c549-979a-05cd-034c-03c2ecc5e3b6"
      },
      "source": [
        "Ok, I surprised myself. This worked out pretty well. I think, I can build an extremely efficient feature from this.\n",
        "\n",
        "But this method is NOT robust yet.\n",
        "\n",
        "- It is not finding the tips, but the points with the greatest distance from center. (look at leaf#19)\n",
        "- It will miserably fail on a more complex, or unfortunately rotated leaf. (look at leaf#78)\n",
        "\n",
        "I really have to smooth the contour to find the dominant symmetry. Also, it would give me a better generalization about the shape. Then I could measure the distance between the smooth contour and the actual one to quantify the variability of the contour.\n",
        "\n",
        "There is a feature in `scipy.ndimage` called 'Mathematical Morphology' that can do stuff like this: \n",
        "\n",
        "![Scipy's Mathematical Morphology functions][1]\n",
        "\n",
        "\n",
        "  [1]: http://www.scipy-lectures.org/_images/morpho_mat.png"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "41ab6867-8a57-a651-7bf6-966879bec502"
      },
      "outputs": [],
      "source": [
        "def cont(img):\n",
        "    return max(measure.find_contours(img, .8), key=len)\n",
        "\n",
        "# let us set the 'brush' to a 6x6 circle\n",
        "struct = [[ 0., 0., 1., 1., 0., 0.],\n",
        "          [ 0., 1., 1., 1., 1., 0.],  \n",
        "          [ 1., 1., 1., 1., 1., 1.], \n",
        "          [ 1., 1., 1., 1., 1., 1.], \n",
        "          [ 1., 1., 1., 1., 1., 1.], \n",
        "          [ 0., 1., 1., 1., 1., 0.],\n",
        "          [ 0., 0., 1., 1., 0., 0.]]\n",
        "\n",
        "erosion = cont(ndi.morphology.binary_erosion(img, structure=struct).astype(img.dtype))\n",
        "closing = cont(ndi.morphology.binary_closing(img, structure=struct).astype(img.dtype))\n",
        "opening = cont(ndi.morphology.binary_opening(img, structure=struct).astype(img.dtype))\n",
        "dilation = cont(ndi.morphology.binary_dilation(img, structure=struct).astype(img.dtype))\n",
        "\n",
        "plt.imshow(img.T, cmap='Greys', alpha=.2)\n",
        "plt.plot(erosion[::,0], erosion[::,1], c='b')\n",
        "plt.plot(opening[::,0], opening[::,1], c='g')\n",
        "plt.plot(closing[::,0], closing[::,1], c='r')\n",
        "plt.plot(dilation[::,0], dilation[::,1], c='k')\n",
        "#plt.xlim([220, 420])\n",
        "#plt.ylim([250, 420])\n",
        "plt.xlim([0, 400])\n",
        "plt.ylim([400, 800])\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "29f4ab22-d829-05c2-7069-a56992ebbb35"
      },
      "source": [
        "This does not make ANY sense. It looks like there is noise around the edges of the original image, and I was not aware of that."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "1ee579b7-d1cb-cfb8-b537-38fbf8da9653"
      },
      "outputs": [],
      "source": [
        "plt.imshow(img.astype(bool).astype(float), cmap='hot')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "a83a7092-cc32-b2ca-46d8-901232c80874"
      },
      "source": [
        "Here we go. We have noise. This challenge just got way more realistic. (Sh*t.) It seems I must learn how to denoise.\n",
        "\n",
        "Ok. The pixels are 8bit greyscale in the original image. Check if I can separate the leaf from the noise by their value.\n",
        "\n",
        "If not, I have a nasty plan. I invert the whole image, do a binary_closing on it to fill up the small holes, then invert it back - now without the noise."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "061a19e1-b1f1-34f7-2aad-32fc0ca89367"
      },
      "outputs": [],
      "source": [
        "erosion = cont(ndi.morphology.binary_erosion(img > 254, structure=struct).astype(img.dtype))\n",
        "closing = cont(ndi.morphology.binary_closing(img > 254, structure=struct).astype(img.dtype))\n",
        "opening = cont(ndi.morphology.binary_opening(img > 254, structure=struct).astype(img.dtype))\n",
        "dilation = cont(ndi.morphology.binary_dilation(img > 254, structure=struct).astype(img.dtype))\n",
        "\n",
        "plt.imshow(img.T, cmap='Greys', alpha=.2)\n",
        "plt.plot(erosion[::,0], erosion[::,1], c='b')\n",
        "plt.plot(opening[::,0], opening[::,1], c='g')\n",
        "plt.plot(closing[::,0], closing[::,1], c='r')\n",
        "plt.plot(dilation[::,0], dilation[::,1], c='k')\n",
        "plt.xlim([0, 400])\n",
        "plt.ylim([400, 800])\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "5a705d21-7134-4596-7042-f09915e99cd6"
      },
      "source": [
        "From these 2 binary morphology tests, it is clear:\n",
        "\n",
        "- the leaf has debris around it's edge\n",
        "- there are not 100% white pixels at the edge\n",
        "\n",
        "But for me this second red contour is just perfect. I will use that as baseline. (yes, I loose some minimal contour information, but it makes my measurements more accurate in exchange)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "792f199c-7694-6bf5-1f42-b334311d5294"
      },
      "outputs": [],
      "source": [
        "# While I was looking up how to threshold an image, I have found this.\n",
        "# The further away from the edges a pixel is, the higher value it gets.\n",
        "# This is important, because it describes the morphology of the leaf better, \n",
        "# than a simple euclidean distance from the center, because it considers\n",
        "# concave parts differently, and that's an important feature I wish to keep.\n",
        "\n",
        "# This is very promising. Using this, I will probably be able to find symmetry.\n",
        "\n",
        "# I also believe, that using this distance map, I will be able to separate\n",
        "# the core shape of the leaf and the edge texture, which are two distinct,\n",
        "# pretty good features.\n",
        "\n",
        "dist_2d = ndi.distance_transform_edt(img)\n",
        "plt.imshow(img, cmap='Greys', alpha=.2)\n",
        "plt.imshow(dist_2d, cmap='plasma', alpha=.2)\n",
        "plt.contour(dist_2d, cmap='plasma')\n",
        "plt.show()"
      ]
    }
  ],
  "metadata": {
    "_change_revision": 0,
    "_is_fork": false,
    "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.0"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}