Analysis/Fit_Ellipse/notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What is the image moment (or how to fit and ellipse)\n",
    "\n",
    "A moment (M) is a specific quantitative measure of a given distribution. The distribution can be one or multi-dimensional ([wikipedia definition](https://en.wikipedia.org/wiki/Moment_(mathematics)).\n",
    "\n",
    "Now let's try to have an intuitive view of what it is in the case of an image.\n",
    "\n",
    "Remember that an image is only a 2D distribution of a set of points. Each point being localized by two coordinates $x$ and $y$. So fo each pixel we have the following $f(x, y) = I_{x, y}$. Where $I_{x, y}$ corresponds to the intensity at the coordinate $x$ and $y$. \n",
    "\n",
    "Before going to intuition let's do some easy math. The formal definition of the image moment (2D and continuous) is:\n",
    "\n",
    "$$M_{pq} = \\int_{-\\infty}^{\\infty} x^p y^q f(x, y) dx dy$$\n",
    "\n",
    "And the 2D discrete equation is : \n",
    "\n",
    "$$M_{pq} = \\sum_x \\sum_y x^p y^q I(x,y)$$\n",
    "\n",
    "So we see here the moment (M) depends on two parameters $p$ and $q$. They are called **the order** of the moment. Different order will give a different **measure** of the distribution we are looking at.\n",
    "\n",
    "For example let's define a simple 1D distribution of length $N$ $f(x) = A$.\n",
    "\n",
    "- the zero order ($p=0$) moment of $A$ is **the number of element** : $M_0 = \\sum_x x^0 = len(A)$ (for pythonist)\n",
    "- the first order ($p=1$) moment of a $A$ if the **sum** of all the element : $M_1 = \\sum_x x^1 = np.sum(A)$ (for pythonist)\n",
    "\n",
    "Now if we divide the first order moment by the zero order moment we have the mean : \n",
    "\n",
    "$$\\bar{A} = \\frac{1}{N} \\sum_x x = \\frac{M_1}{M_0}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The image moment \n",
    "\n",
    "Now let's try to apply this 1D example to an image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Do some import\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from tifffile import TiffFile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQwAAAEACAYAAABGTkjoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEbBJREFUeJzt3W+sHNV9xvHvs9cYMFD+lMqG4MhWhCOI2gJJCEkgKS2h\nECW0rwpIiRCgvKENRFUJNi/6NggpSvqGF00BUQQ0DgkIK1HABqSmispfOxhfHEITGtxgA0mgNJFa\n8P76Ymbt9Y3vvWfvnZmdM/N8JMt75+7ds39mnz3nzNn5KSIwM0sxmPYdMLN8ODDMLJkDw8ySOTDM\nLJkDw8ySOTDMLNmigSFpk6RdknZKulfSkZJOkrRV0ouSHpF0wpzr/0TSbkkX13v3zaxJWmgdhqR1\nwGPAGRHxv5K+CXwP+ADwRkTcKukm4MSI2CjpTOBe4MPAe4BtwIaIGNb7MMysCYv1MP4beAdYJWkF\nsAr4BXAZcFd5nbuAvywv/wVwX0S8ExEvAy8B51Z9p81sOhYMjIj4FfBV4OcUQfFmRGwFVkfEvvJq\n+4DV5eVTgT1jN7GHoqdhZh2wYGBIeh/wJWAdRRgcK+lz49eJYkyz0Ppyrz0364gVi/z+Q8API+KX\nAJK+A3wU2CtpTUTslXQK8Fp5/f8C1o79/WnltkNIcoiYTUlEaKl/u9gcxm7gPElHSxJwETALbAGu\nKq9zFfBgefkh4ApJKyWtB04HnlzqnTOzdlmwhxERP5L0z8DTwBB4FvhH4Dhgs6RrgZeBvyqvPytp\nM0WovAtcF/46rFlnLHhYtbZGPSQxm5o6hyRmZgc4MMwsmQPDzJI5MMwsmQPDzJI5MMwsmQPDzJI5\nMMwsmQPDzJI5MMwsmQPDzJI5MMwsmQPDzJI5MMwsmQPDzJI5MMwsmQPDzJI5MMwsmQPDzJKl1FZ9\nv6TtY//eknS966ua9c9EJwGWNKCoM3Iu8EWWWF/VJwE2m54mTwJ8EfBSRLyC66ua9c6kgXEFcF95\n2fVVzXomOTAkrQQ+C3xr7u9cX9WsHybpYVwKPBMRr5c/75O0BmAp9VXNLD+TBMaVHByOQFFH1fVV\nzXok6SiJpGOA/wTWR8Tb5baTgM3Aeynrq0bEm+XvbgauoaivekNEPDzn9jxEMZuS5RwlcW1Vs55x\nbVUza4QDw8ySOTDMLJkDw8ySOTDMLJkDw8ySOTDMLJkDw8ySOTDMLJkDw8ySOTDMLJkDw8ySOTDM\nLJkDw8ySOTDMLJkDw8ySOTDMLJkDw8ySOTDMLFlSYEg6QdL9kl6QNCvpI66tatY/qT2MfwC+FxFn\nAH8E7AY2AlsjYgPwaPkzZW3Vy4EzgUuA28qarGaWuZTq7ccDF0TEHQAR8W5EvIVrq5r1Tson/3rg\ndUl3SnpW0jfKOiWurWrWMymBsQI4B7gtIs4BfkM5/BhxbVWzfkgJjD3Anoh4qvz5fooA2evaqmb9\nsmhgRMRe4BVJG8pNFwG7gC24tqpZr6xIvN4XgXskrQT+A7gamAE2S7qWsrYqQETMStoMzFLUVr0u\nplGP0cwq59qqZj3j2qpm1ggHhpklc2CYWTIHhpklc2CYWTIHhpklc2CYWbLUhVudNhwOD7s9IpiZ\nmWn43pi1Vy8DYxQQ0sLrVyQxvrBtOBw6QKzXehMYVaxoHQwGh72diGAw8OjOuq+zgTEcDhftQVRl\nbk9k1L57I9Y1nfxYjIjGwmI+o97IfPMjZjnqTA+jyR7FJOb2Ptp4H81SZd/DGL0hc3kjRoR7Hpat\nrHsYkrJ9483teXji1HKQbWC0dQiyVOMB0qXHZd2S7Udal99Uo2GLWdtkGRhdDotxDg5rmywDY//+\n/dO+C43yJKm1RWpt1ZclPSdpu6Qny21Tq63alx7GuNEch4PDpim1hxHAn0TE2RExKnvo2qpTkPOR\nIcvfJG/kuR/rrq06Je5t2LRM0sPYJulpSV8ot7m26pQ5OKxpqeswPh4Rr0r6A2CrpN3jv4yIWKTW\niKf6azQeHP7Cm9UpqYcREa+W/78OPEAxxNjn2qrtMvrCW9+OIllzFg0MSaskHVdePga4GNhJUUPV\ntVVbaL7zdpgtV8qQZDXwQHkocwVwT0Q8IulpXFu11bzU3KqWZW1VH1pcOn/JzXpXW9UdlqXzkRVb\njmy/rWrL4xP72FJk2cOw6o2+6OYjLLaQbAPD4/B6+FykthC/6+ywPNdhh5PtHIYnPpvhuQ4bl3UP\nwztv83wS437LOjBsejxk6ScHhi3LKDg8ROyH7APD385sDwdH92UfGO4St4/XdHRX9oFh7eVvzXZP\nJwLDw5J2c2h0RycCw8OS9nNodEMnAsPy4NDInwPDGuXQyFtnAsNfRsuHQyNfnXmXeSfMi1+vPHUm\nMCw/Do38pNZWnSnrqm4pf55aXdVF7meTzVkFHBp5Se1h3EBxFvDRq+u6qlYZh0Y+UuqSnAZ8Gvgn\nDtZXbW1dVfcy8uRl5HlI+fT/GnAjML46qtV1Vb3yMz9eRp6HBQND0meA1yJiO79bvR0o6qqycO3U\nxvcCr/zMl1+7dlvsFH0fAy6T9GngKOD3JN1NWVc1Iva2ta7q3FPLWR4ksX//fvcSW2rBHkZE3BwR\nayNiPXAF8FhEfJ5M6qp6p8uTF+G116QnAR59ZN9CBnVV3b01q1aWtVUnMRgMPAOfKR/xqkfvaqtO\nwr2MfPm1a5/OBwb4kypXft3apxeBAcXO50nQ/Hg42S69CQxwFzdHPmLSLr17NdzNNVu63gUGODRy\n4wV47dHLwLD8ODTaobeB4V6G2eR6GxiWH/cypq/XgeFeRn58pGu6eh0Y4NDIjV+v6ep9YIB3wtx4\nMdf0ODBKDo18eDHX9PiZH+PQyId7GdPhwJjDn1558GkLpsPvjjkiwl9Sy4TDvXl+xg/Dh+7y4bUZ\nzXJgzMPzGfnw0KQ5DowFODTy4PmM5ixWl+QoSU9I2iFpVtJXyu2trK1aB4dGHhwazVj0JMCSVkXE\nbyWtAP4N+DuKUolvRMStkm4CToyIjWVt1XuBD1NUPNsGbIiI4ZzbzGrg6Z0xHw74xdV6EuCI+G15\ncSUwA/yaFtdWrcNwOPSMfCY8CVqvlGLMA0k7KGqoPh4Ru2h5bdU6eEfMh3uD9UnpYQwj4iyKsoef\nkHThnN+3rrZqXdzdzYOHkPVJ7mdHxFvAd4EPUtZWBWhrbdW6ODTy4CFkPRY7SnLy6AiIpKOBTwHb\nyaS2al28EjQPHkZWb7HaqqcAd0kaUITL3RHxqKTtZFBbtS5eCZoPV4KvVudrq9apg1nYSR5GHsq1\nVafEn1x58ARodRwYyzAcDh0aGfAEaHX8TC7TcDhEknfKlvPwsRreyysSEQ6OlnNoLJ/37oqNgsOs\nixwYNZHk+Y0Wci9jeRwYNfKX1qxrvDfXzEOU9nEvY+kcGA3xhKh1gffgBkWEQ8Oy5r23YR6itIOH\nJUvjwJgSH0WxHDkwpmi0StSmw986npwDowUcGtPh531yDoyW8M47He5lTMaB0SKSHBwN8/M9GQdG\nC3kntrZyYLSUQ6M5PsSazoHRYg6N5visXGlSChmtlfS4pF2Snpd0fbm9N/VVp8nzGs3wCtw0KbVV\n1wBrImKHpGOBZyhKI17NEuurduUkwE1zgZ569SWY666tujcidpSX/wd4gSIIelVftQ18DtF6+RDr\n4ibqh0laB5wNPEEP66u2gVeH1sfP6+KSA6McjnwbuCEi3h7/XZ/qq7aFd26bhqTAkHQERVjcHRGj\nsoi9ra/aFj7HRvV8iHVhKUdJBNwOzEbE18d+1ev6qm3hc2xYk1KOkpwP/CvwHAeHFpsoQmAz8F7K\n+qoR8Wb5NzcD11DUV70hIh6ec5uO8Rr407EaXR/uLecoiWurdoxDoxpdDg3XVrUDfGIeq5MDo4N8\n6NXq4sDoMIeGVc2B0XEeoliVHBg94CGKVcWB0SMODVsuB0bPeHWoLYf3nB4aFVNycNikvMf0mIPD\nJuU9xQ4Eh4+mFLxadn4ODDtgdDSl7+HhHtf8/MzYYfX17F4+krQwB4bNq2+h4aHI4hwYtqC+hIbP\nK5LGz5Atquuh4bBI52fJknR1ebnDYjJ+pmwiXVq34bCYnJ8tm1gX1m04LJYm5STAd0jaJ2nn2DaX\nSbQDw5Tc3niDwSC7+9wWKc/ancAlc7ZtBLZGxAbg0fJnyjKJlwNnln9zmyS/Mh036nHkYDAY+PDp\nMqSUSvwB8Os5m10m0X7HaJVoWz+9HRbLt9RX1mUSbV6jHkebwiMiHBYVWPar6TKJtpDx8JjWJKkn\nOKuz1GfRZRJtYqNJ0iaDwxOc1VrqM+kyibZkTQWH5yyqt2KxK0i6D/gkcLKkV4C/B24BNku6lrJM\nIkBEzEraDMxSlEm8LvyK2TzGV49KYjgcVnbbDot6uFSitcpgMGD//v0Hfp6ZmTnk5xS5HOKdFtdW\ntc4ZvelH++fcIDkcT26mcWBYb8wXHO5VpHMxZuuN4XB4SC9iMBg4LBq06KSnWdvktBS9a9zDMLNk\nDgwzS+bAMLNkDgwzS+bAMLNkDgwzS+bAMLNkDgwzS+bAMLNkDgwzS+bAMLNkDgwzS+bAMLNkDgwz\nS+bAMLNktQSGpEvK2qo/kXRTHW2YWfMqP0WfpBngx8BFFDVJngKujIgXxq7jU/SZTUnbTtF3LvBS\nRLwcEe8A/0JRc9XMMldHYLwHeGXsZ9dXNeuIOs7pmTrceLaGtudzCvCq23N7bo9zlvPHdQTG3Pqq\nazm0ovvIsu74Epzi9tye21ueOiY9V1BMev4Z8AuK2qqHTHqaWZ4q72FExLuS/gZ4GJgBbndYmHXD\nVCqfmVmevNLTzJI1Hhh1rAKVdIekfZJ2jm07SdJWSS9KekTSCWO/21S2v1vSxRO2tVbS45J2SXpe\n0vU1t3eUpCck7ZA0K+krdbY3dhszkrZL2lJ3e5JelvRc2d6TDbR3gqT7Jb1QPqcfqfH1e3/5uEb/\n3pJ0fc2Pb1O5f+6UdK+kIytrLyIa+0cxp/ESsA44AtgBnFHB7V4AnA3sHNt2K/Dl8vJNwC3l5TPL\ndo8o78dLwGCCttYAZ5WXj6WY4D2jrvbK21hV/r8C+Hfg/DrbK2/nb4F7gIfqfD7L2/gZcNKcbXW2\ndxdwzdhzenzdz2d5WwOKw6dra9w/1wE/BY4sf/4mcFVV7TUdGB8Fvj/280ZgY0W3vY5DA2M3sLq8\nvAbYXV7eBNw0dr3vA+cto90HKZbB194esIpiqf0H6mwPOA3YBlwIbKn7+aQIjN+fs62W9ijC4aeH\n2d7E63cx8IOaH99JFB9iJ1KE4RbgU1W11/SQpMlVoKsjYl95eR+wurx8KoeuC1nyfZC0jqJn80Sd\n7UkaSNpR3u7jEbGrzvaArwE3AsOxbXW2F8A2SU9L+kLN7a0HXpd0p6RnJX1D0jE1tjfuCuC+8nIt\n7UXEr4CvAj+nWNbwZkRsraq9pgNjKodkoojOhdqe+H5JOhb4NnBDRLxdZ3sRMYyIsyg++T8h6cK6\n2pP0GeC1iNgOHPZLSjU8nx+PiLOBS4G/lnRBje2toFg0eFtEnAP8hqKnW1d7AEhaCXwW+Nbv3Fi1\nr9/7gC9R9LhPBY6V9Lmq2ms6MFJXgVZhn6Q1AJJOAV6b5z6cVm5LJukIirC4OyIerLu9kYh4C/gu\n8MEa2/sYcJmkn1F8Gv6ppLtrbI+IeLX8/3XgAYovMNbV3h5gT0Q8Vf58P0WA7K359bsUeKZ8jFDf\n4/sQ8MOI+GVEvAt8h2IqoJLH13RgPA2cLmldmbiXAw/V1NZDFJM9lP8/OLb9CkkrJa0HTqdYjZpE\nkoDbgdmI+HoD7Z08mtGWdDTFeHR7Xe1FxM0RsTYi1lN0oR+LiM/X+PhWSTquvHwMxTh/Z42Pby/w\niqQN5aaLgF0UY/3K2xtzJQeHI6PbraO93cB5ko4u99WLgFmqenxLmbxZzj+KpP0xxWzspopu8z6K\n8dr/UcyRXE0x+bMNeBF4BDhh7Po3l+3vBv58wrbOpxjb76B4424HLqmxvT+k+KLeDuA54MY4OLlV\neXtz2v4kB4+S1PX41pePbQfw/GifqPPxAX9MMXn8I4pP4ONrbu8Y4A3guLFtdbb3ZYoQ3ElxROiI\nqtrzSk8zS+aVnmaWzIFhZskcGGaWzIFhZskcGGaWzIFhZskcGGaW7P8BaFChZWtPn/4AAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f30b03ddf60>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Load the image\n",
    "tf = TiffFile(\"binary_cell.tif\")\n",
    "\n",
    "# Get the numpy array\n",
    "a = tf.asarray()\n",
    "\n",
    "# Replace all 255 to 1 so the image is now made of \"0\" and \"1\"\n",
    "a[a == 255] = 1\n",
    "print(np.unique(a))\n",
    "\n",
    "_ = plt.imshow(a, interpolation=\"none\", cmap=\"gray\", origin='lower')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal here is to use the image moment to **fit an ellipse**. The ellipse can be defined by four parameters : the center, the orientation, the major axis and the minor axis.\n",
    "\n",
    "Note it exists many different method to fit an ellipse.\n",
    "\n",
    "### Find the centroid\n",
    "\n",
    "First we are going to detect the centroid (the center of the ellipse) of our object by computing the zero and first order of the image.\n",
    "\n",
    "Remember the discrete equation : \n",
    "\n",
    "$$M_{pq} = \\sum_x \\sum_y x^p y^q I(x,y)$$\n",
    "\n",
    "Since we are working on a binary image (only $0$ and $1$) values we remove $I(x,y)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The zero order moment is then :\n",
    "\n",
    "$$M_{00} = \\sum_x \\sum_y x^0 y^0$$\n",
    "\n",
    "\n",
    "Let's do it in python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "142369"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The sum of all value of 1\n",
    "M_00 = np.sum(a)\n",
    "M_00"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's compute the first order moment for $x$ :\n",
    "\n",
    "$$M_{10} = \\sum_x \\sum_y x^1 y^0$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "46616805"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Here we get all the coordinates of pixel equal to 1\n",
    "xx, yy = np.where(a == 1)\n",
    "\n",
    "M_10 = np.sum(xx)\n",
    "M_10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's compute the first order moment for $y$ :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "48021475"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M_01 = np.sum(yy)\n",
    "M_01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the centroid $C$ is given by :\n",
    "\n",
    "\n",
    "$$C_x = \\bar{x} = \\frac{M_{10}}{M_{00}}$$\n",
    "\n",
    "$$C_y = \\bar{y} = \\frac{M_{01}}{M_{00}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C = (327.43648547085388, 337.30288897161603)\n"
     ]
    }
   ],
   "source": [
    "C_x = M_10 / M_00\n",
    "C_y = M_01 / M_00\n",
    "print(\"C =\", (C_x, C_y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's verify it visually :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQwAAAEACAYAAABGTkjoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEhFJREFUeJzt3V2sXNV5xvH/M8cYMFCMS2VDcGorwhFEbYEkhCSQlNah\nECWkVwWkRAhQbmgDUVWCzUXDXRFSlPSGi6aAKAKKQwLFShRsA2oTReXTDsYHB2hCgxtsIAmUJlIL\nnrcXe489Pvics2Zm7z374/lJlufsM2fWfOx5Zq211+xXEYGZWYretO+AmTWHA8PMkjkwzCyZA8PM\nkjkwzCyZA8PMki0aGJI2StolaaekuyUdKWmFpK2Snpe0RdLyOdd/QdJuSReUe/fNrEpaaB2GpDXA\nI8BpEfG/ku4Fvgd8AHg9Im6WdD1wQkRskHQ6cDfwYeA9wDZgXUT0y30YZlaFxXoY/w28DSyTtARY\nBvwCuBi4I7/OHcCf55c/B9wTEW9HxEvAi8DZRd9pM5uOBQMjIn4FfA34OVlQvBERW4GVEbEvv9o+\nYGV++WRgz9BN7CHraZhZCywYGJLeB3wZWEMWBsdK+vzwdSIb0yy0vtxrz81aYskiv/8Q8KOI+CWA\npO8AHwX2SloVEXslnQS8ml//v4DVQ39/Sr7tEJIcImZTEhEa928Xm8PYDZwj6WhJAtYDs8Bm4PL8\nOpcDD+SXHwQulbRU0lrgVODxce+cmdXLgj2MiPixpH8CngT6wNPAPwDHAZskXQW8BPxFfv1ZSZvI\nQuUd4Orw12HNWmPBw6qlNeohidnUlDkkMTM7wIFhZskcGGaWzIFhZskcGGaWzIFhZskcGGaWzIFh\nZskcGGaWzIFhZskcGGaWzIFhZskcGGaWzIFhZskcGGaWzIFhZskcGGaWzIFhZskcGGaWLKW26vsl\nbR/696aka1xf1ax7RjoJsKQeWZ2Rs4EvMWZ9VZ8E2Gx6qjwJ8HrgxYh4GddXNeucUQPjUuCe/LLr\nq5p1THJgSFoKfBb41tzfub6qWTeM0sO4CHgqIl7Lf94naRXAOPVVzax5RgmMyzg4HIGsjqrrq5p1\nSNJREknHAP8JrI2It/JtK4BNwHvJ66tGxBv5724AriSrr3ptRDw05/Y8RDGbkkmOkri2qlnHuLaq\nmVXCgWFmyRwYZpbMgWFmyRwYZpbMgWFmyRwYZpbMgWFmyRwYZpbMgWFmyRwYZpbMgWFmyRwYZpbM\ngWFmyRwYZpbMgWFmyRwYZpbMgWFmyRwYZpYsKTAkLZd0n6TnJM1K+ohrq5p1T2oP4++B70XEacAf\nAruBDcDWiFgHPJz/TF5b9RLgdOBC4Ja8JquZNVxK9fbjgfMi4jaAiHgnIt7EtVXNOiflk38t8Jqk\n2yU9LembeZ0S11Y165iUwFgCnAXcEhFnAb8hH34MuLaqWTekBMYeYE9EPJH/fB9ZgOx1bVWzblk0\nMCJiL/CypHX5pvXALmAzrq1q1ilLEq/3JeAuSUuB/wCuAGaATZKuIq+tChARs5I2AbNktVWvjmnU\nYzSzwrm2qlnHuLaqmVXCgWFmyRwYZpbMgWFmyRwYZpbMgWFmyRwYZpYsdeFWq/X7/cNujwhmZmYq\nvjdm9dXJwBgEhLTw+hVJDC9s6/f7DhDrtM4ERhErWnu93mFvJyLo9Ty6s/ZrbWD0+/1FexBFmdsT\nGbTv3oi1TSs/FiOisrCYz6A3Mt/8iFkTtaaHUWWPYhRzex91vI9mqRrfwxi8IZvyRowI9zyssRrd\nw5DU2Dfe3J6HJ06tCRobGHUdgoxrOEDa9LisXRr7kdbmN9Vg2GJWN40MjDaHxTAHh9VNIwNj//79\n074LlfIkqdVFam3VlyQ9I2m7pMfzbVOrrdqVHsawwRyHg8OmKbWHEcAfR8SZETEoe+jaqlPQ5CND\n1nyjvJHnfqy7tuqUuLdh0zJKD2ObpCclfTHf5tqqU+bgsKqlrsP4eES8Iun3gK2Sdg//MiJikVoj\nnuov0XBw+AtvVqakHkZEvJL//xpwP9kQY59rq9bL4AtvXTuKZNVZNDAkLZN0XH75GOACYCdZDVXX\nVq2h+c7bYTaplCHJSuD+/FDmEuCuiNgi6UlcW7XWvNTcitbI2qo+tDg+f8nNOldb1R2W8fnIik2i\nsd9Wtcn4xD42jkb2MKx4gy+6+QiLLaSxgeFxeDl8LlJbiN91dlie67DDaewchic+q+G5DhvW6B6G\nd97q+STG3dbowLDp8ZClmxwYNpFBcHiI2A2NDwx/O7M+HBzt1/jAcJe4frymo70aHxhWX/7WbPu0\nIjA8LKk3h0Z7tCIwPCypP4dGO7QiMKwZHBrN58CwSjk0mq01geEvozWHQ6O5WvMu807YLH69mqk1\ngWHN49BontTaqjN5XdXN+c9Tq6u6yP2ssjkrgEOjWVJ7GNeSnQV88Oq6rqoVxqHRHCl1SU4BPg38\nIwfrq9a2rqp7Gc3kZeTNkPLp/3XgOmB4dVSt66p65WfzeBl5MywYGJI+A7waEdt5d/V2IKurysK1\nUyvfC7zys7n82tXbYqfo+xhwsaRPA0cBvyPpTvK6qhGxt651VeeeWs6aQRL79+93L7GmFuxhRMQN\nEbE6ItYClwKPRMQXaEhdVe90zeRFePU16kmABx/ZN9GAuqru3poVq5G1VUfR6/U8A99QPuJVjs7V\nVh2FexnN5deuflofGOBPqqby61Y/nQgMyHY+T4I2j4eT9dKZwAB3cZvIR0zqpXOvhru5ZuPrXGCA\nQ6NpvACvPjoZGNY8Do166GxguJdhNrrOBkan3HjjtO9BIdzLmL7Wr/RcTCd2Qgla8jgjwkdOJuSV\nnhPw0KRZ/HpNV+cDA7wTNo0Xc02PAyPn0GgOD0mmx8/8kFaExo03ZnMWw//g3dsaPhHqXsZ0dH7S\ncy5J7VtC3qJJz2H9ft/fDxqDJz0LFBHeCRvCQ5Pq+Rk/jNb1MFqsE4fFa8SBMY9WzGd0hOczquPA\nWEBrQuOrX532PSiVT8NYnQUnPSUdBfwrcCSwFPiXiNgoaQVwL/D75CcBjog38r/ZCFwJ7AeuiYgt\nh7ndRvUj3e1tBk+Cpplk0nPRoySSlkXEbyUtAX4I/A1ZqcTXI+JmSdcDJ0TEhry26t3Ah8kqnm0D\n1kVEf85tNuod6E+w5mhNr7BEpR4liYjf5heXAjPAr6lxbdUy9Pt9z8g3hHuD5UopxtyTtIOshuqj\nEbGLmtdWLYN3xOZwb7A8KT2MfkScQVb28BOSzp/z+9rVVi2Lu7vN4CFkeZL72RHxJvBd4IPktVUB\n6lpbtSwOjWbwELIci1VvP1HS8vzy0cCngO00pLZqWTwT3wweRhZvsdqqJwF3SOqRhcudEfGwpO00\noLZqWbwStDlcCb5Y/vLZBFqYha3kYeSh/OWzKfEnVzN4ArQ4DowJeGVhM3gCtDh+JifU7/eR5J2y\n5jx8LIb38oJEhIOj5hwak/PeXbBBcJi1kQOjJJI8v1FD7mVMxoFRIn9pzdrGe3PJPESpH/cyxufA\nqIgnRK0NvAdXyHVBrem891bMQ5R68LBkPA6MKfFRFGsiB8YUDVaJ2nT4W8ejc2DUgENjOvy8j86B\nURPeeafDvYzRODBqRJKDo2J+vkfjwKgh78RWVw6MmnJoVMeHWNM5MGrMoVEdn5UrTUoho9WSHpW0\nS9Kzkq7Jt6+QtFXS85K2DM4unv9uo6QXJO2WdEGZD6DtPK9RDa/ATZNSW3UVsCoidkg6FniKrDTi\nFYxZX7UtJwGumgv0lKsrwVx2bdW9EbEjv/w/wHNkQdCp+qp14HOIlsuHWBc3Uj9M0hrgTOAxOlhf\ntQ68OrQ8fl4XlxwY+XDk28C1EfHW8O+6VF+1Lrxz2zQkBYakI8jC4s6IGJRF7Gx91brwOTaK50Os\nC0s5SiLgVmA2Ir4x9KtO11etC59jw6qUcpTkXODfgGc4OLTYSBYCm4D3ktdXjYg38r+5AbiSrL7q\ntRHx0JzbdIyXwJ+OxWj7cG+SoySurdoyDo1itDk0XFvVDvCJeaxMDowW8qFXK4sDo8UcGlY0B0bL\neYhiRXJgdICHKFYUB0aHODRsUg6MjvHqUJuE95wOGhRTcnDYqLzHdJiDw0blPcUOBIePpmS8WnZ+\nDgw7YHA0pevh4R7X/PzM2GF19exePpK0MAeGzatroeGhyOIcGLagroSGzyuSxs+QLartoeGwSOdn\nyZK0dXm5w2I0fqZsJG1at+GwGJ2fLRtZG9ZtOCzGk3IS4Nsk7ZO0c2ibyyTagWFK0954vV6vcfe5\nLlKetduBC+ds2wBsjYh1wMP5z+RlEi8BTs//5hZJfmVabtDjaIJer+fDpxNIKZX4A+DXcza7TKK9\ny2CVaF0/vR0Wkxv3lXWZRJvXoMdRp/CICIdFASZ+NV0m0RYyHB7TmiT1BGdxxn0WXSbRRjaYJK0y\nODzBWaxxn0mXSbSxVRUcnrMo3pLFriDpHuCTwImSXgb+FrgJ2CTpKvIyiQARMStpEzBLVibx6vAr\nZvMYXj0qiX6/X9htOyzK4VKJViu9Xo/9+/cf+HlmZuaQn1M05RDvtLi2qrXO4E0/2D/nBsnheHIz\njQPDOmO+4HCvIp2LMVtn9Pv9Q3oRvV7PYVGhRSc9zeqmSUvR28Y9DDNL5sAws2QODDNL5sAws2QO\nDDNL5sAws2QODDNL5sAws2QODDNL5sAws2QODDNL5sAws2QODDNL5sAws2QODDNLVkpgSLowr636\ngqTry2jDzKpX+Cn6JM0APwHWk9UkeQK4LCKeG7qOT9FnNiV1O0Xf2cCLEfFSRLwN/DNZzVUza7gy\nAuM9wMtDP7u+qllLlHFOz9ThxtMltD2fk4BX3J7bc3ucNckflxEYc+urrubQiu4DE93xMZzk9tye\n25tMGZOeS8gmPf8U+AVZbdVDJj3NrJkK72FExDuS/gp4CJgBbnVYmLXDVCqfmVkzeaWnmSWrPDDK\nWAUq6TZJ+yTtHNq2QtJWSc9L2iJp+dDvNubt75Z0wYhtrZb0qKRdkp6VdE3J7R0l6TFJOyTNSvq7\nMtsbuo0ZSdslbS67PUkvSXomb+/xCtpbLuk+Sc/lz+lHSnz93p8/rsG/NyVdU/Lj25jvnzsl3S3p\nyMLai4jK/pHNabwIrAGOAHYApxVwu+cBZwI7h7bdDHwlv3w9cFN++fS83SPy+/Ei0BuhrVXAGfnl\nY8kmeE8rq738Npbl/y8B/h04t8z28tv5a+Au4MEyn8/8Nn4GrJizrcz27gCuHHpOjy/7+cxvq0d2\n+HR1ifvnGuCnwJH5z/cClxfVXtWB8VHg+0M/bwA2FHTbazg0MHYDK/PLq4Dd+eWNwPVD1/s+cM4E\n7T5Atgy+9PaAZWRL7T9QZnvAKcA24Hxgc9nPJ1lg/O6cbaW0RxYOPz3M9ipevwuAH5T8+FaQfYid\nQBaGm4FPFdVe1UOSKleBroyIffnlfcDK/PLJHLouZOz7IGkNWc/msTLbk9STtCO/3UcjYleZ7QFf\nB64D+kPbymwvgG2SnpT0xZLbWwu8Jul2SU9L+qakY0psb9ilwD355VLai4hfAV8Dfk62rOGNiNha\nVHtVB8ZUDslEFp0LtT3y/ZJ0LPBt4NqIeKvM9iKiHxFnkH3yf0LS+WW1J+kzwKsRsR047JeUSng+\nPx4RZwIXAX8p6bwS21tCtmjwlog4C/gNWU+3rPYAkLQU+CzwrXfdWLGv3/uAL5P1uE8GjpX0+aLa\nqzowUleBFmGfpFUAkk4CXp3nPpySb0sm6QiysLgzIh4ou72BiHgT+C7wwRLb+xhwsaSfkX0a/omk\nO0tsj4h4Jf//NeB+si8wltXeHmBPRDyR/3wfWYDsLfn1uwh4Kn+MUN7j+xDwo4j4ZUS8A3yHbCqg\nkMdXdWA8CZwqaU2euJcAD5bU1oNkkz3k/z8wtP1SSUslrQVOJVuNmkSSgFuB2Yj4RgXtnTiY0ZZ0\nNNl4dHtZ7UXEDRGxOiLWknWhH4mIL5T4+JZJOi6/fAzZOH9niY9vL/CypHX5pvXALrKxfuHtDbmM\ng8ORwe2W0d5u4BxJR+f76npglqIe3ziTN5P8I0van5DNxm4s6DbvIRuv/R/ZHMkVZJM/24DngS3A\n8qHr35C3vxv4sxHbOpdsbL+D7I27HbiwxPb+gOyLejuAZ4Dr4uDkVuHtzWn7kxw8SlLW41ubP7Yd\nwLODfaLMxwf8Ednk8Y/JPoGPL7m9Y4DXgeOGtpXZ3lfIQnAn2RGhI4pqzys9zSyZV3qaWTIHhpkl\nc2CYWTIHhpklc2CYWTIHhpklc2CYWbL/B7zL0FwmKkHiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f30b0537390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(C_x, C_y, color='red', marker='+', s=100)\n",
    "_ = plt.imshow(a, interpolation=\"none\", cmap=\"gray\", origin='lower')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well it seems to be exact !\n",
    "\n",
    "### Find the major and minor axis\n",
    "\n",
    "Here it becomes a little bit tricky.\n",
    "\n",
    "The information about the object orientation can be derived using the **second order central moments** to construct a **covariance matrix** (see [Wikipedia](https://en.wikipedia.org/wiki/Image_moment) for more details).\n",
    "\n",
    "Here is a new concept : **central moment** ($\\mu$) which describes **distribution of mean** (unlike moment (M) which describes distribution only).\n",
    "\n",
    "Note that sometime moment is also called **raw moment**.\n",
    "\n",
    "The discrete equation of the central moment is the : \n",
    "\n",
    "$$\\mu_{pq} = \\sum_{x} \\sum_{y} (x - \\bar{x})^p(y - \\bar{y})^q f(x,y)$$\n",
    "\n",
    "with $\\bar{x}=\\frac{M_{10}}{M_{00}}$ and $\\bar{y}=\\frac{M_{01}}{M_{00}}$\n",
    "\n",
    "Now it becomes difficult to get the intuition of what's going next. From [Wikipedia](https://en.wikipedia.org/wiki/Image_moment) :  \n",
    "\n",
    "> The covariance matrix of the image $I(x,y)$ is :\n",
    "\n",
    "> $$\\operatorname{cov}[I(x,y)] = \\begin{bmatrix} \\mu'_{20}  & \\mu'_{11} \\\\ \\mu'_{11} & \\mu'_{02} \\end{bmatrix}$$\n",
    "\n",
    "> The eigenvectors of the covariance matrix of the image correspond to the **major** and **minor** axes of the image intensity, so the orientation can thus be extracted from the **angle of the eigenvector** associated with the **largest eigenvalue**. It can be shown that this angle Θ is given by the following formula:\n",
    "\n",
    "> $$\\Theta = \\frac{1}{2} \\arctan \\left( \\frac{2\\mu'_{11}}{\\mu'_{20} - \\mu'_{02}} \\right)$$\n",
    "\n",
    "> Where :\n",
    "\n",
    "> $$\\mu'_{20} = \\mu_{20} / \\mu_{00} = M_{20}/M_{00} - \\bar{x}^2$$\n",
    "> $$\\mu'_{02} = \\mu_{02} / \\mu_{00} = M_{02}/M_{00} - \\bar{y}^2$$\n",
    "> $$\\mu'_{11} = \\mu_{11} / \\mu_{00} = M_{11}/M_{00} - \\bar{x}\\bar{y}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first compute the second order raw moment ($M_{20}$, $M_{02}$ and $M_{11}$) : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "M_20 = np.sum(xx ** 2)\n",
    "M_02 = np.sum(yy ** 2)\n",
    "M_11 = np.sum(xx * yy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute $\\mu'_{20}$, $\\mu'_{02}$ and $\\mu'_{11}$ :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu_20 = M_20 / M_00 - C_x ** 2\n",
    "mu_02 = M_20 / M_00 - C_y ** 2\n",
    "mu_11 = M_11 / M_00 - C_x * C_y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the **orientation** $\\Theta$ : \n",
    "\n",
    "$$\\Theta = \\frac{1}{2} \\arctan \\left( \\frac{2\\mu'_{11}}{\\mu'_{20} - \\mu'_{02}} \\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "angle = -36.72080335329064°\n"
     ]
    }
   ],
   "source": [
    "theta = 1/2 * np.arctan((2 * mu_11) / (mu_20 - mu_02))\n",
    "\n",
    "# Convert it in degree\n",
    "angle = np.rad2deg(theta)\n",
    "print(\"angle = {}°\".format(angle))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the eigenvalues with this equation : \n",
    "\n",
    "$$\\Delta = \\sqrt{4{\\mu'}_{11}^2 + ({\\mu'}_{20}-{\\mu'}_{02})^2}$$\n",
    "$$\\lambda_i = \\frac{\\mu'_{20} + \\mu'_{02}}{2}  \\pm \\frac{\\Delta}{2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "delta = np.sqrt(4 * mu_11 ** 2 + (mu_20 - mu_02) ** 2)\n",
    "\n",
    "lambda_1 = ((mu_20 + mu_02) + delta) / 2\n",
    "lambda_2 = ((mu_20 + mu_02) - delta) / 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the major and minor axes from the eigenvalues : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Major axis = 590.4638690697747\n",
      "Minor axis = 139.87004447009642\n"
     ]
    }
   ],
   "source": [
    "semi_major_length = np.sqrt(np.abs(lambda_1)) * 2\n",
    "semi_minor_length = np.sqrt(np.abs(lambda_2)) * 2\n",
    "\n",
    "major_length = semi_major_length * 2\n",
    "minor_length = semi_minor_length * 2\n",
    "\n",
    "print(\"Major axis = {}\".format(major_length))\n",
    "print(\"Minor axis = {}\".format(minor_length))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively we can also compute (with Numpy) the **orientation** and **axes length** from the eigenvectors of the covariance matrix.\n",
    "\n",
    "$$\\operatorname{cov}[I(x,y)] = \\begin{bmatrix} \\mu'_{20}  & \\mu'_{11} \\\\ \\mu'_{11} & \\mu'_{02} \\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cov = np.asarray([[mu_20, mu_11],\n",
    "                  [mu_11, mu_02]])\n",
    "eigvalues, eigvectors = np.linalg.eig(cov)\n",
    "\n",
    "# Get the associated eigenvectors and eigenvalues\n",
    "eigval_1, eigval_2 = eigvalues\n",
    "eigvec_1, eigvec_2 = eigvectors[:, 0], eigvectors[:, 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the orientation from the **first eigenvector** with $ \\Theta = \\operatorname{atan2}(y, x) \\quad$ and the **axes length** from the **eigenvalues**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "angle = -36.72080335329064°\n",
      "Major axis = 590.4638690697747\n",
      "Minor axis = 139.87004447009653\n"
     ]
    }
   ],
   "source": [
    "theta = np.arctan2(eigvec_1[1], eigvec_1[0])\n",
    "angle = np.rad2deg(theta)\n",
    "print(\"angle = {}°\".format(angle))\n",
    "\n",
    "semi_major_length = np.sqrt(np.abs(eigval_1)) * 2\n",
    "semi_minor_length = np.sqrt(np.abs(eigval_2)) * 2\n",
    "\n",
    "major_length = semi_major_length * 2\n",
    "minor_length = semi_minor_length * 2\n",
    "print(\"Major axis = {}\".format(major_length))\n",
    "print(\"Minor axis = {}\".format(minor_length))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We find the **same value** as above. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's see how our vector looks (which is supposed to define the major orientation of our object)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQwAAAEACAYAAABGTkjoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFHhJREFUeJzt3X+wXGV9x/H3Z3MTIAm/Ip0EJDYZJQqObUAEVIQKiGDV\nOmNHcJRxhHHGQQR/VEmiU7B2wHGK2qHSsVaooEQjSITRCuFHLY5VCCQQSMKPShSEBER+CRaS7Ld/\nnLNhc0nuffbePXv2nPN5zWSye+69+5z99dnnec6z56uIwMwsRavsHTCz6nBgmFkyB4aZJXNgmFky\nB4aZJXNgmFmycQND0mJJd0taI+lySbtImiVphaR7JV0naa9Rv3+fpPWSji92981skDTWOgxJ84Ab\ngQMj4nlJ3wd+ArwW+H1EfFnS2cDeEbFI0kHA5cAbgJcD1wMLIqJd7N0ws0EYr4fxNLAZmC5pBJgO\nPAy8G/h2/jvfBt6TX/4bYGlEbI6IDcD9wGH93mkzK8eYgRERfwAuAH5LFhRPRsQKYHZEbMp/bRMw\nO7+8H/BQ1008RNbTMLMaGDMwJL0S+AQwjywMZkr6YPfvRDamGWt9udeem9XEyDg/PxT4RUQ8DiDp\nh8AbgY2S5kTERkn7Ao/mv/87YG7X3++fb9uOJIeIWUkiQhP92/HmMNYDR0jaTZKA44C1wDXAh/Lf\n+RCwPL98NXCypGmS5gMHALdMdOfMbLiM2cOIiDskXQqsBNrA7cC/AbsDyySdBmwA3pf//lpJy8hC\nZQtwevjrsGa1MeZh1cIa9ZDErDRFDknMzLZxYJhZMgeGmSVzYJhZMgeGmSVzYJhZMgeGmSVzYJhZ\nMgeGmSVzYJhZMgeGmSVzYJhZMgeGmSVzYJhZMgeGmSVzYJhZMgeGmSVzYJhZMgeGmSVLqa36akmr\nuv49JelM11c1a56eTgIsqUVWZ+Qw4ONMsL6qTwJsVp5BngT4OOD+iHgQ11c1a5xeA+NkYGl+2fVV\nzRomOTAkTQPeBfxg9M9cX9WsGXrpYZwI3BYRj+XXN0maAzCR+qpmVj29BMb7eXE4AlkdVddXNWuQ\npKMkkmYAvwHmR8Qz+bZZwDLgFeT1VSPiyfxnS4BTyeqrnhUR1466PQ9RzEoymaMkrq1q1jCurWpm\nA+HAMLNkDgwzS+bAMLNkDgwzS+bAMLNkDgwzS+bAMLNkDgwzS+bAMLNkDgwzS+bAMLNkDgwzS+bA\nMLNkDgwzS+bAMLNkDgwzS+bAMLNkDgwzS5YUGJL2knSFpHWS1ko63LVVzZontYfxz8BPIuJA4C+A\n9cAiYEVELABuyK+T11Y9CTgIOAG4KK/JamYVl1K9fU/gLRFxMUBEbImIp3BtVbPGSfnknw88JukS\nSbdL+mZep8S1Vc0aJiUwRoBDgIsi4hDgWfLhR4drq5o1Q0pgPAQ8FBG35tevIAuQja6tatYs4wZG\nRGwEHpS0IN90HHA3cA2urWrWKCOJv/dx4LuSpgH/C3wYmAIsk3QaeW1VgIhYK2kZsJasturpUUY9\nRjPrO9dWNWsY11Y1s4FwYJhZMgeGmSVzYJhZMgeGmSVzYJhZMgeGmSVLXbhVa+12e4fbI4IpU6YM\neG/MhlcjA6MTENLY61ck0b2wrd1uO0Cs0RoTGP1Y0dpqtXZ4OxFBq+XRndVfbQOj3W6P24Pol9E9\nkU777o1Y3dTyYzEiBhYWO9PpjexsfsSsimrTwxhkj6IXo3sfw7iPZqkq38PovCGr8kaMCPc8rLIq\n3cOQVNk33uiehydOrQoqGxjDOgSZqO4AqdP9snqp7Edand9UnWGL2bCpZGDUOSy6OThs2FQyMLZu\n3Vr2LgyUJ0ltWKTWVt0g6U5JqyTdkm8rrbZqU3oY3TpzHA4OK1NqDyOAv4qIgyOiU/bQtVVLUOUj\nQ1Z9vbyRR3+su7ZqSdzbsLL00sO4XtJKSR/Jt7m2askcHDZoqesw3hwRj0j6M2CFpPXdP4yIGKfW\niKf6C9QdHP7CmxUpqYcREY/k/z8GXEU2xNjk2qrDpfOFt6YdRbLBGTcwJE2XtHt+eQZwPLCGrIaq\na6sOoZ2dt8NsslKGJLOBq/JDmSPAdyPiOkkrcW3Voeal5tZvlayt6kOLE+cvuVnjaqu6wzJxPrJi\nk1HZb6va5PjEPjYRlexhWP91vujmIyw2lsoGhsfhxfC5SG0sftfZDnmuw3aksnMYnvgcDM91WLdK\n9zD84h08n8S42SodGFYeD1mayYFhk9IJDg8Rm6HygeFvZw4PB0f9VT4w3CUePl7TUV+VDwwbXv7W\nbP3UIjA8LBluDo36qEVgeFgy/Bwa9VCLwLBqcGhUnwPDBsqhUW21CQx/Ga06HBrVVZt3mV+E1eLn\nq5pqExhWPQ6N6kmtrTolr6t6TX69tLqq4+znIJuzPnBoVEtqD+MssrOAd55d11W1vnFoVEdKXZL9\ngXcA/86L9VWHtq6qexnV5GXk1ZDy6f9V4DNA9+qooa6r6pWf1eNl5NUwZmBIeifwaESs4qXV24Gs\nripj104d+KvAKz+ry8/dcBvvFH1vAt4t6R3ArsAeki4jr6saERuHta7q6FPLWTVIYuvWre4lDqkx\nexgRsSQi5kbEfOBk4MaIOIWK1FX1i66avAhvePV6EuDOR/aXqEBdVXdvzfqrkrVVe9FqtTwDX1E+\n4lWMxtVW7YV7GdXl52741D4wwJ9UVeXnbfg0IjAge/F5ErR6PJwcLo0JDHAXt4p8xGS4NO7ZcDfX\nbOIaFxjg0KgaL8AbHo0MDKseh8ZwaGxguJdh1rvGBkajnHtu2XvQF+5llK/2Kz3H04gXoQQ1uZ8R\n4SMnk+SVnpPgoUm1+PkqV+MDA+r9Inxh6wtl70LfeTFXeRwYuTqGRjvaLLhwAa8+o15DLw9JyuNH\nvkstQuPcc7M5C4lWawrf+JffcO8+0PqHFg/srW0/q/pEqHsZ5Wj8pOdokmq3hPzpXcWei7PLF554\nIWccdka5O9Qn7Xbb3w+agMlMejowdqB259DIj5KcctUpfOfO7zB3j7ls+MQGWjWoAFGLXuGAOTAK\nUKcxf/dh1RsfuJFjLz0WgPs+fh+vmvWqMvesLxwavfFh1QLU9UV4zPxjeGbxMwAccOEBfOV/vlLy\nHk1erXqDQ86BMYbahMY552x3dea0mcQ5wWkHn8anr/s0s/9pNlvb1X3T1W4IOcTGHJJI2hX4GbAL\nMA34UUQsljQL+D7w5+QnAY6IJ/O/WQycCmwFzoyI63Zwu5Xq79dqeDLKzb+5maP+4ygA1n1sHa/Z\n5zUl79HEeRI0TaFzGJKmR8RzkkaAnwN/R1Yq8fcR8WVJZwN7R8SivLbq5cAbyCqeXQ8siIj2qNus\n1Duw7p9gz21+jhnnzQDg/GPPZ9GRi0reo4mrTa+wQIXOYUTEc/nFacAU4AmGuLZqEdrtdq0XC02f\nOp04Jzj90NNZfMNi9jh/D7a0t5S9WxNS597gMEgpxtyStJqshupNEXE3Q15btQhNeCF+/a+/zi9P\n+yXPvPAMU784lbsevavsXZqQOvcGy5bSw2hHxEKysodHSXrrqJ8PXW3VojShu3v4/ofzp8/9CYDX\n/evr+MJ/faHkPepd3YeQZUruZ0fEU8CPgdeT11YFGNbaqkVpQmjsOrIrcU7wqSM+xbk/O5epX5zK\n5q2by96tntR5CFmm8aq37yNpr/zybsDbgFVUpLZqUZoyE3/B2y/g1o/cypb2Fqb94zTu2HhH2bvU\nkyYMIwdtvBjeF7gxn8P4FXBNRNxAVlv1bZLuBY7JrxMRa4FObdX/pOTaqkWp23dNxnLofofy/Oef\nB2DhNxay5IYl2/38vJvP4+nnny5j15J4aNJfXho+CTXMwjF97obPcd7PzwPg+c8/z5XrruQDV36A\nk157Ekv/dmnJe7dzTRhG9sLfJSlJEyfX7th4Bwu/sXC7bdOnTufS91zKew96b0l7NTYv6NqeA6NE\nTQyNP77wR3Y/f/fttu2xyx7cc8Y9zJk5p6S9Gpt7GS/yl89K1G63kdSoWflPXvtJdhvZbbttz21+\njpOvOHloh2nDul9V05xXecEiohHBsfLhlVyy6pKXbN/S3sLKh1dy0cqLStirNA6NyfOQpCB1fnE+\n/tzjLF+/nItXX8xtD9/GSGuEZzc/C8D0kems+ugqFrxsQcl7uWMemngOY2g1YX7jiT89wfL1y7lk\n9SXc8rtb2NzezIH7HMjqj65mpDVS9u7tUNNDw4ExxOp4jtCdefL/nuRH63/E9+76HkfPO3pov/Xq\nwHBgDL06D1F25NkXnmXGtBll78ZONTk0HBgV0aTexrBzYExMvaf0h4zrglrV+dU7YJ3Dr1aupg0R\n+8WBURJJXq5slePAKFFnlaiVw/NJvXNgDAGHRjn8uPfOgTEk/OIth3sZvXFgDBFJDo4B8+PdGwfG\nEPKL2IaVA2NIOTQGx4dY0zkwhphDY3Dq/iXBfkkpZDRX0k2S7pZ0l6Qz8+2zJK2QdK+k6zpnF89/\ntljSfZLWSzq+yDtQd57XGAyvwE2TUlt1DjAnIlZLmgncRlYa8cNMsL5qU79LMllN+Lp8mZoSzEXX\nVt0YEavzy38E1pEFQaPqqw4Dn8y2WD7EOr6e+mGS5gEHk9UoaVx91WHg1aHF8eM6vuTAyIcjVwJn\nRcQz3T9rUn3VYeEXt5UhKTAkTSULi8siolMWsbH1VYdFE046PGg+xDq2lKMkAr4FrI2Ir3X9qNH1\nVYeFz7Fhg5RylORI4L+BO3lxaLGYLASWAa8ANgDvi4gn879ZApwKbCEbwlw76jYd4wXwp2N/1H24\n51P02TYOjf6oc2j4FH22jU/MY0VyYNSQD71aURwYNebQsH5zYNSchyjWTw6MBvAQxfrFgdEgDg2b\nLAdGw3h1qE2GXzkN1Cmm5OCwXvkV02AODuuVXym2LTh8NCXj1bI758CwbTpHU5oeHu5x7ZwfGduh\npp7dy0eSxubAsJ1qWmh4KDI+B4aNqSmh4fOKpPEjZOOqe2g4LNL5UbIkdV1e7rDojR8p60md1m04\nLHrnR8t6Vod1Gw6LiUk5CfDFkjZJWtO1zWUSbdswpWpvvFarVbl9HhYpj9olwAmjti0CVkTEAuCG\n/Dp5mcSTgIPyv7lIkp+Zmuv0OKqg1Wr58OkkpJRKvBl4YtRml0m0l+isEh3WT2+HxeRN9Jl1mUTb\nqU6PY5jCIyIcFn0w6WfTZRJtLN3hUdYkqSc4+2eij6LLJFrPOpOkgwwOT3D210QfSZdJtAkbVHB4\nzqL/Rsb7BUlLgaOBfSQ9CPw98CVgmaTTyMskAkTEWknLgLVkZRJPDz9jthPdq0cl0W63+3bbDoti\nuFSiDZVWq8XWrVu3XZ8yZcp211NU5RBvWVxb1Wqn86bvvD5HB8mOeHIzjQPDGmNnweFeRToXY7bG\naLfb2/UiWq2Ww2KAxp30NBs2VVqKXjfuYZhZMgeGmSVzYJhZMgeGmSVzYJhZMgeGmSVzYJhZMgeG\nmSVzYJhZMgeGmSVzYJhZMgeGmSVzYJhZMgeGmSVzYJhZskICQ9IJeW3V+ySdXUQbZjZ4fT9Fn6Qp\nwD3AcWQ1SW4F3h8R67p+x6foMyvJsJ2i7zDg/ojYEBGbge+R1Vw1s4orIjBeDjzYdd31Vc1qoohz\neqYON24voO2d2Rd4xO25PbfHIZP54yICY3R91blsX9G9Y1I7PgH7uj235/Ymp4hJzxGySc9jgYfJ\naqtuN+lpZtXU9x5GRGyRdAZwLTAF+JbDwqweSql8ZmbV5JWeZpZs4IFRxCpQSRdL2iRpTde2WZJW\nSLpX0nWS9ur62eK8/fWSju+xrbmSbpJ0t6S7JJ1ZcHu7SvqVpNWS1ko6v8j2um5jiqRVkq4puj1J\nGyTdmbd3ywDa20vSFZLW5Y/p4QU+f6/O71fn31OSziz4/i3OX59rJF0uaZe+tRcRA/tHNqdxPzAP\nmAqsBg7sw+2+BTgYWNO17cvAZ/PLZwNfyi8flLc7Nd+P+4FWD23NARbml2eSTfAeWFR7+W1Mz/8f\nAX4JHFlke/ntfAr4LnB1kY9nfhsPALNGbSuyvW8Dp3Y9pnsW/Xjmt9UiO3w6t8DX5zzg18Au+fXv\nAx/qV3uDDow3Aj/tur4IWNSn257H9oGxHpidX54DrM8vLwbO7vq9nwJHTKLd5WTL4AtvD5hOttT+\ntUW2B+wPXA+8Fbim6MeTLDBeNmpbIe2RhcOvd7B9EM/f8cDNBd+/WWQfYnuTheE1wNv61d6ghySD\nXAU6OyI25Zc3AbPzy/ux/bqQCe+DpHlkPZtfFdmepJak1fnt3hQRdxfZHvBV4DNAu2tbke0FcL2k\nlZI+UnB784HHJF0i6XZJ35Q0o8D2up0MLM0vF9JeRPwBuAD4LdmyhicjYkW/2ht0YJRySCay6Byr\n7Z73S9JM4ErgrIh4psj2IqIdEQvJPvmPkvTWotqT9E7g0YhYBezwS0oFPJ5vjoiDgROBj0l6S4Ht\njZAtGrwoIg4BniXr6RbVHgCSpgHvAn7wkhvr7/P3SuATZD3u/YCZkj7Yr/YGHRipq0D7YZOkOQCS\n9gUe3ck+7J9vSyZpKllYXBYRy4turyMingJ+DLy+wPbeBLxb0gNkn4bHSLqswPaIiEfy/x8DriL7\nAmNR7T0EPBQRt+bXryALkI0FP38nArfl9xGKu3+HAr+IiMcjYgvwQ7KpgL7cv0EHxkrgAEnz8sQ9\nCbi6oLauJpvsIf9/edf2kyVNkzQfOIBsNWoSSQK+BayNiK8NoL19OjPaknYjG4+uKqq9iFgSEXMj\nYj5ZF/rGiDilwPs3XdLu+eUZZOP8NQXev43Ag5IW5JuOA+4mG+v3vb0u7+fF4Ujndotobz1whKTd\n8tfqccBa+nX/JjJ5M5l/ZEl7D9ls7OI+3eZSsvHaC2RzJB8mm/y5HrgXuA7Yq+v3l+Ttrwfe3mNb\nR5KN7VeTvXFXAScU2N7ryL6otxq4E/hMvDi51ff2RrV9NC8eJSnq/s3P79tq4K7Oa6LI+wf8Jdnk\n8R1kn8B7FtzeDOD3wO5d24ps77NkIbiG7IjQ1H6155WeZpbMKz3NLJkDw8ySOTDMLJkDw8ySOTDM\nLJkDw8ySOTDMLNn/A/n31udvLf2LAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f30b0268358>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "\n",
    "# Plot the centroid in red\n",
    "plt.scatter(C_x, C_y, color='red', marker='+', s=100)\n",
    "\n",
    "# Plot the first eigenvector\n",
    "scale = 100\n",
    "x1, x2 = [C_x, eigvec_1[0] * scale]\n",
    "y1, y2 = [C_y, eigvec_1[1] * scale]\n",
    "plt.arrow(x1, y1, x2, y2, color='green', lw=1, head_width=20)\n",
    "\n",
    "# Show the image\n",
    "_ = plt.imshow(a, interpolation=\"none\", cmap=\"gray\", origin='lower')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **orientation** seems to be correct."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we draw the major and minor axes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQwAAAEACAYAAABGTkjoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXd4VVXWxt8VwIIoCgkhSBWBAREFBLsMiogN7GDFMuB8\nyiifDbGO3zgzTiiiMlgoiggoIE1BqSJKL6GGiFRpCb1DSHLe7499IiGmnHvvKfcm6/c868m9J+fs\ndco+7919CUkoiqI4IS7oE1AUJXZQwVAUxTEqGIqiOEYFQ1EUx6hgKIriGBUMRVEcU6xgiEhPEVkt\nIitFZISInC4ilURkmoisFZGpInJuvv1/FZE0EWnr7ekriuInUtQ4DBGpDWAmgIYkM0XkKwCTAVwE\nYDfJZBHpAeA8ki+LSCMAIwC0AHA+gOkA6pO0vL0MRVH8oLgSxkEAWQDKi0hZAOUBbAfQHsBQe5+h\nAO6wP3cAMJJkFslNANYBaOn2SSuKEgxFCgbJvQD6APgNRij2k5wGIJFkhr1bBoBE+3M1AFvzJLEV\npqShKEoJoEjBEJG6ALoDqA0jBhVE5KG8+9DUaYoaX65jzxWlhFC2mP9fBmAuyT0AICJjAVwJIF1E\nqpJMF5EkADvt/bcBqJHn+Or2tlMQERURRQkIkhLuscW1YaQBuEJEzhQRAdAGQCqAbwB0tvfpDGC8\n/XkigE4icpqI1AFQD8DCcE9OUZToosgSBsnlIvI5gMUALABLAXwC4GwAo0TkCQCbANxn758qIqNg\nRCUbwFPU6bCKUmIoslvVM6daJVGUwPCySqIoivI7KhiKojhGBUNRFMeoYCiK4hgVDEVRHKOCoSiK\nY1QwFEVxjAqGoiiOUcFQFMUxKhiKojhGBUNRFMeoYCiK4hgVDEVRHKOCoSiKY1QwFEVxjAqGoiiO\nUcFQFMUxKhiKojhGBUNRFMc4ia3aQERS8tgBEXlG46sqSukjpEWARSQOJs5ISwB/Q5jxVXURYEUJ\nDj8XAW4DYB3JLdD4qopS6ghVMDoBGGl/1viqilLKcCwYInIagNsBjM7/P42vqiilg1BKGDcDWEJy\nl/09Q0SqAkA48VUVRYk9QhGM+3GyOgKYOKoaX1VRShGOeklE5CwAmwHUIXnI3lYJwCgANWHHVyW5\n3/7fKwAeh4mv+izJKfnS0yqKogREJL0kGltVUUoZGltVURRfUMFQFMUxKhiKojhGBUNRFMeoYCiK\n4hgVDEVRHKOCoSiKY1QwFEVxjAqGoiiOUcFQFMUxKhiKojhGBUNRFMeoYCiK4hgVDEVRHKOCoSiK\nY1QwFEVxjAqGoiiOUcFQFMUxKhiKojjGkWCIyLkiMkZE1ohIqohcrrFVFaX04bSE8R6AySQbAmgC\nIA3AywCmkawPYIb9HXZs1Y4AGgFoB2CAHZNVUZQYx0n09ooAriU5BABIZpM8AI2tqiilDie//HUA\n7BKRT0VkqYgMtOOUaGxVRSllOBGMsgCaARhAshmAI7CrH7lobFVFKR04EYytALaSXGR/HwMjIOka\nW1VRShfFCgbJdABbRKS+vakNgNUAvoHGVlWUUkVZh/v9DcBwETkNwHoAjwEoA2CUiDwBO7YqAJBM\nFZFRAFJhYqs+xSDiMSqK4joaW1VRShkaW1VRFF9QwVAUxTEqGIqiOEYFQ1EUx6hgKIriGBUMRVEc\no4KhKIpjnA7cKtFYllXgdpIoU6aMz2ejKNFLqRSMXIEQKXr8iogg78A2y7JUQJRSTakRDDdGtMbF\nxRWYDknExWntTin5lFjBsCyr2BKEW+QvieT619KIUtIokT+LJH0Ti8LILY0U1j6iKLFIiSlh+Fmi\nCIX8pY9oPEdFcUrMlzByX8hYeRFJaslDiVliuoQhIjH74uUveWjDqRILxKxgRGsVJFzyCkhJui6l\nZBGzghHSS7VxIzBnDpCWBuzbB2RlAWeeCVSpAlSrBjRqZOzss7074RBQ4VCilZhccctxVWTBAuDF\nF41QtG4N/OlPQEICUK4ccOwYkJEBbNsGpKYCa9YA1asDf/6zsVatjJhEASociptEsuJWTAqGk+rI\nlk96o/rrvSDJycCDDwJliylM5eQAq1YBs2YZmz0bqFEDuPNO4K67gMaNgQBfXG3jUNwiEsH4vdW+\nKINZ5HcFgBQAC+1tlQBMA7AWwFQA5+bZvyeAX2FCKrYtID1GYsWRvX4d951Vli/2vokHjh8odv+C\nE8kmZ88mu3cna9YkL7yQfO01ct268NJzCcuyIrp3amp08M4XZo5KGCKyEUBzknvzbEsGsJtksoj0\nAHAeyZft2KojALSAiXg2HUB9klaeY4t3WgTFnvMLLyBLgG7XHcLMTTPx1T1foVlSs0gcAkuWAF98\nAYwYYao2jz0G3HsvUKFC+OmGfTpa2lDChz6UMDYCqJxvWxpMuEQAqAogjSdLFz3y7Pc9gCvyHRup\nQhZN8+bknDkkyZErRzI+OZ4fLPiAlmWF86N+KpmZ5NixZPv25Hnnkd26kWlpkacbBlraUAvHGEEJ\nw+nPFAFMF5HFItLF3ha9sVV37DDtDwA6Ne6EuY/PxZCUIbh39L3Yf3x/ZGmfdppp15gwAVixAqhY\nEbjuOqBdO2DyZMDHcSG5XbGxOhZFiT2cCsbVJJsCuBnA0yJybd5/klEWW/W004ATJ37/Wq9yPcx9\nYi6qVqiKZh83w6Jti4o4OASqVwfefhvYvBm4/37gjTeAhg2BwYNP8e81ucKRk5Pjm0+ldOJIMEju\nsP/uAjAOQEsAGVEbW7VWLWD9+lM2nVH2DPS/pT+Sb0zGrSNuxXvz33NlyrtJ/Aygc2dg0SLg44+B\nUaOAunWBd98FDh92x4cDcie8qXAoXlGsYIhIeRE52/58FoC2AFbCxFDtbO/WGdEUW7VlS2BhwS7v\naXQP5j0xD8NWDMOdX92Jvcf2FrhfWIiYMRxTpgDjxwNz5wJ16gD/+Adw8KB7foqhsHU7FCViimvk\nAFAHwDLbVgHoyZPdqtNRcLfqKwDWwTSM3lRAmpE22hTN5MnkVVcVucvxrON89rtnWevdWpy3ZV7x\naYZLWhr50ENkQgL5n/+QR45456sQIr3faiXLGEGjZ9gHRuQ0wgsWkaLfkMxMsnJlcvPmYl+mcWvG\nMSE5gb3m9GKOlVPs/mGzahV5zz1kUhL5/vvk8ePe+SoC7VlRo9fjMNzGjWDMxZ73U0+ZYeBvvVVs\nWpv2b0KnMZ0QXz4eQ+8YisrlK0d6eoWzdKlpHF2xAnj9dTOeo7hRqB5AHctRaqHX4zDcNrijkkWz\nZo2pBhw96uiXNzM7k89PeZ41+tbgz5t/dnRMRMydS7ZqRTZqZKpQbowRiQA3nolabBhLW5XEvuji\nuf12U/wPgYlpE1mlVxX++6d/e1tFIY1ITJhA1q9P3ngjuXy5t/4ckJOTE3iGVvPWWBoFo9h2DJJc\nupRMTCQPHgzppdm8fzOvGnwV233RjjsP7wzp2LA4cYL84AOyShXyL38hd+zw3mcxaFtHyTVG8O6W\n7Eps06bADTcAffqEdFjNijUxq/MsXJJ4CZp+3BSzN8/26ARtypUDunUDfvkFOPdc4KKLgH/+Ezh+\n3Fu/RaCjSJUCiURtwjW4p5TFs3Gj6TFZvz7E31jD5LWTmdgrkf/48R/MzskOK42QWb+e7NDBzJCd\nPNkfnw5x69mpBWeM4N2N2V6SXByd/zvvmPUtJk0Ka02LbQe34f6v78fpZU/HF3d+gcQKicUf5AaT\nJwPPPANcfLEZNVq7tj9+HUBqL0uswgh6SUrHE3/uOTPfY9SosA4//5zzMbPzTFx+/uVo9kkzzNw4\n0+UTLIRbbjGL+jRrBjRvbuatBFhNyYtWWUopkRRPwjW4W7xyxoIFplFx2zbnxxTA1HVTmdQ7iW/+\n8KZ/VRSS3LDhZDVlyhT//IaAm89VTaskv+NmlSQuLs75ZKu33jLzO777DoigOL3j0A48OPZBEMSI\nu0Yg6eyksNMKmcmTgaefBq69FujbF4iP98+3Q3QN0uiGpblKElKR+JVXgP37gQEDIvKZdHYSpj08\nDa1qtUKzT5ph2vppEaUXErfcAqxcaYSicWNg+HAgANEvitxfI501WwKJpHgSrsH9IpZz1q4l4+PN\nGA0XmLFhBqv1qcZXZ7zKrJwsV9J0zMKFZJMmZLt25KZN/voOAbeft1pwVZKYL2EACC1Ker16QP/+\nwD33mBglEXJ9neuxtOtSLNi2ANcPvR7bDvq49EeLFsDixWbFr+bNgX79zOrnUQajrASkhE/Mt2Hk\nEvJ1dO9uFtmZMCGi9oxccqwcvPPzO/hg4Qf47I7P0O7CdhGnGRJr1wJduwKZmcDQoUD9+v76d4C2\nbUQHLI2Tz/JbyGRmmjUz3n479GOL4MdNP7J63+rsMa0HT2SfcDXtYsnJIfv3N1Wud98136MML569\nmn9VktJbwgBM1LOWLYEPPwTat3ftXHYd2YVHxj+Cg5kHMfLukahZsaZraTti3Trg0UeBMmWATz8F\nLrjAX//FoCWNYKGWMBxORiuIBQvML3JKSnjHF0KOlcN3fnqHVXpV4cS0ia6m7YjsbLJ3b3NtH34Y\n+PT5/HiRB9S8L2GUGMGwb0R4fPUVWaNGxIO6CuLnzT+zRt8afO7755iZnel6+sWSmkq2aGGmz//2\nm//+iyDoF6e0GiN4d0tEL0nE3HefaTDs0AE4etTVpK+ueTVSnkzB2r1rcd2n12HT/k2upl8sDRua\nwWqtWpmelNGj/fVfBAygOqxEiBNVAVAGJq7qN/b3sOOq0sMSBiIpZVgW+cgj5G23kVnuj6ewLIu9\n5/RmQnICx60Z53r6jli40Awtf+wx8tChYM6hALzMD2rBlDCeBZBqOwSAlwFMI1kfwAz7O+y4qh0B\nNALQDsAAEYmNUowIMGgQkJ1tShsu//qJCJ6/6nlMvH8iun/fHc9+9ywyszNd9VEsLVoAKSnmWps2\nNXFUogBqSSN2KE5RYAIRTQfQGidLGGHHVaXHJQxEUsogzS9vy5bkyy9Hlk4R7D26lx1GdmDzj5tz\n/d7w1umImFGjzGS8f/3LNJAGjC4NWHJKGO8CeBFA3kkb0RtXFSGO/MxPhQpm3Yxx48zISQ8478zz\nMK7jODzc5GFcMegKjEkd44mfIrn3XjNKdMoUsyrZli3+n0MeNPhSbFCkYIjIbQB2kkwBUGDfLfl7\nHbQwfM8FEa/REB9vXqS+fc04Bg8QETx7xbOY9MAkvDTtJTw96Wkcz/Z5rYsaNYAZM4CbbjLVlcmT\n/fVfALq+RpRTVPEDwL8AbAGwEcAOAEcADIOpklS190nCySrJywBe5qlVkssLSNevoldkpKWR1aqR\nX3wReVpFsO/YPt791d1s+lFTrt291lNfhfLTT8yqXo2zO11lFiUOEK2eRG+VJJSxE61wsg0jGXZb\nhS0S79ifG8GEVDwNJsTietjzVRiAYMTFxUWee1etIqtWNXV+D7Esi/0X9Gd8cjxHrhzpqa+CyMrJ\n4u3vXcFfW9Yjr76a3LLF93PIS9AvVUk2+igYE+3PYcdVpY+CYd+cyFm2zDQQjh/vTnpFsGT7El74\n/oXsOrErj55wFoTJDV6b8RrbfN6G2VknTENoYiL53Xe++c9P0C9VSTb6IRhump83x5VSBkkuWmQi\nqU2a5E56RXDg+AF2HN2RTT5swrRdaZ77m7Z+Gqv1qcb0Q+knN/74I1m9OvnKK56MS3FC0C9WSTWq\nYBR7g9xh3jzfRMOyLH68+GPGJ8dz2PJhnvnZfnA7k3onccaGGX/8Z0aGGVLeujW504eATvnQYEoq\nGIEIhquiMXeuEY2xY91LswiWpy9ngw8a8PHxj/PIiSOupp2dk83Wn7Xmmz+8WcRO2aaUUbOmKWX5\nTNAvV0k0qmA4M9eqJ0uWmDr+iBHupFcMhzIP8aGxD/Gi/17E1TtXu5bu33/4O1t/1trZ6udjx5qZ\nr4MHu+bfCdpjooIRmGDYN8sdVq4kk5J8e4Esy+LgpYNZ+T+V+WnKpxGnN3PDTCb1TuL2g9udH5Sa\nSjZoQP71r2YBIp8I+gUraUYVjIBEIy3NTIvv39+9NIthZcZKNuzfkI+Me4SHMsObQJZ+KJ3V+lTj\n1HVTQz/4wAETH+XKKz1ZDqAggs4vJc0YwbsbGxPDXMa1FZ8aNAB+/NEEe+7Vy500i6FxlcZY1GUR\n4iQOLQa2wMqMlSEdb9HCw+MexmOXPoYb694Y+gmccw4wdqwJd9CiBTB/fuhphIjRDCUqiERtwjVE\nh8q6x5YtZKNG5PPP+7qO5tBlQxmfHM+BSwbScrii1ts/vs1rh1zrTkiEb74xDcDDh0eelgOCzjMl\nxahVkigQjT17zKLCDz3k69Dq1J2pbDygMR/4+gEePH6wyH1/3PQjE3slcuuBre6dwMqVZO3apifF\nY7EMOr+UFKNWSaKASpWAadNMZLX27YEjR3xx2zChIRb8ZQHKly2PywZehmXpy/6409//jl1HduHB\nsQ/iszs+w/nnuDiBuHFjYMECUzW7915Pr5taNQmeSNQmXEMUqGyuuU5WllnRqmVLctcu99MvguEr\nhjM+OZ4DFg44pYqSI2C7L9qxx7Qe3jk/fpx89FGyaVNP1w7VwVzBljBKvWB4IhqWRfbsabogfQ5h\n+MvuX3jJh5fwvtH3cf+x/STJd64Grxp8lfdxUiyLTE42M3wXLPDMTdD5JdaNKhhRKBok2a+feXl8\nHiF5LOsY//rNX1n3vboctGQQE18Af9vv44rhEyaYxtAJEzxJXgdzqWBEhXnCuHFmhOQ4/xf+HbRk\nEOPeimPX2+C4F8U1Fi40A9sGDPAk+aDzSiwbVTCiXDQWLTIljT59vA8o9OabJEAL4G33g4+3B5t1\nBe+6D9x3BswjB8x+XrNuHVmvnlkf1eUeFC1lqGBEhYUdQa04Nm0iGzc2w6p9mC7ee05vthzYkpnZ\nmTxeBuw2qRvr9KvDBVu9a1sokF27zKjQBx4wDaMuoqKhghEV5toktfzs30+2bUvedJMZYu0R87bM\nY0JyAjfu22g2mEzCMavHMCE5gX3n9vW3inL0KHnnnWaa/L59riYddF6JRaMKhic31RtOnCCffNKU\nNjZscD35vUf3sta7tU4NlpTnejbs3cAWn7Rg+5HtuefoHtf9F0p2NvnMM+TFF5M7driadNB5JdaM\nKhgxJhqWZXpQEhPJGQUsXBN2shY7jOzAZyY/c+o/8l1LZnYmu3/XnbXercW5v811zb+DEyTffpus\nW9dVsdSqiQpG1JinTJ9u1gr94ANXGkP7zevH5h835/GsfG0FhTRwjl8znlV6VWHyz8nMsfybA8P+\n/c3yf6vdW9tDRSMKBAPAGQAWwKwEngrg3/b2qI2t6tEN9o7160315IknImoUXLh1IROSE0KOpLZp\n3yZeMegK3jL8Fu464uPI1GHDTAlr4ULXklTRCFgwaF7u8vbfsgDmA7gGJszAS/b2HvhjmIFyAGrD\nrB4eV0Cagd+0UMyzRtBcDh0yjYJXXRVW/X7fsX2s068OR68eHZb7E9kn+MKUF1ijbw3+tPmnsNII\niwkTzBiVmTNdSzLovBILRj+qJADKA1gE4CJEeWxVL8yz7tZccnLIt94yC/KE8KtrWRbv/upuPvXt\nUxGfwre/fMvEXon81+x/+VdFmTnTiMY337iWZNB5JdqNHpcw4mBKDYcAJNvb9uX5v+R+B/ABgAfz\n/G8QgLsLSDPwmxbmjfaecePMsOqPPnLUrtF/QX9e+tGlPJZ1zBX3Ww5s4dWDr+ZNw25ixuEMV9Is\nlgULXB1KrlWTAAWDJ1/yijBVktbIIxj2//aycMG4q4C0Ar9pEdxs7/nlF9P9+PDD5JHCVwpfsn0J\n45PjXQ+vmJWTxZ7Te/L8Pudz1sZZrqZdKIsWuRosSkUjYMGgedFfB/ACYiS2qoc33HuOHDGC0bix\nEZB8HDh+gBe+f6GnYRW///V7Vu1dlf836/+crSweKYsXG9FwKYRD0PkkWo1eCQaAeNg9IADOBDAb\nwA2IkdiqXpnnjaC5WJapmiQkkF9/nWezxY6jO7LrxK6en8LWA1t53afX8YahN3DHIXcHXBXIkiVG\nNPJcbyQEnVei0eihYFwMYKktAisAvGhvj5nYqh7edP9YtIisVYt87jnyxAl+tOgjNvmwiW+xV7Ny\nsvj6zNdZrU81Tl8/3XuHS5eaLtcxYyJOSqsmPgqGVxb0DYtJ0di9m7z5Zh6+rAmbvnyeLzFX8zNt\n/TQm9U7iGzPf8L6KkpJiShouhKUMOp9Em1EFIxjzrWqSS04O5z97L4+dd44rv77hsOPQDl4/9Hr+\n+bM/c9tBj+OSzJtnulxnzYooGS1lqGBEjfkuGqTphrzgAjNV/qg/1ZK8ZOdk861Zb7Fq76qcsm6K\nt86mTzdtOBGOCA06n0STUQUjePN8YFd+9u8nO3Y0vSguzskIhR82/sBqfaqx5/Se7sQ5KYwJE0yb\nxsqVESUTdB6JFqMKRvSYr8JhWeSgQabYPnCg96t5FUDG4Qy2HdaW1wy5hlsObPHO0fDh5Pnnk5s3\nR5RM0PkjGowqGNFnvrJ6tRnodffdvoc2IMkcK4f/nP1PJvZK5KS1kTdSFkqfPibC3N69YScRdL6I\nBqMKRnSar+0bx46ZUI3VqpHffeef3zzM3jSb1ftW54tTX/QupEH37mSrVhHN7A06XwRtjODdFXP/\n/EVE/HcaECICy7L8c/jDD8CjjwK33moCRJ91ln++Aew+uhuPjHsE+4/vx5f3fImaFWu668CygI4d\ngbJlgeHDgbjQg/e5Fow7RiEZ9g3QUIkeQ9LfDNq6NbB8OXD4MNC0qQlj6CPx5ePx7QPf4o4/3YEW\nA1tg4i8T3XUQFwcMGwZs3Qq89lpYSQTxI1liiKR4Eq4hCoplQZjvPSmjRpnBT2++6WuA6Fzm/DaH\nNd+tyf/9/n+ZmZ3pbuI7d5og0CPDm0sTdF4I0qhVktjB9yrK9u3AE08AGRnAp58Cl1zin28Ae4/t\nxaPjH0XGkQx8efeXqHNeHfcSX74caNMGmDIFaNYspENLc7WEEVRJtIQRnMr7h2WRQ4aYAVBvvEFm\nuvxrX6x7i33n9mVCcgLHprozE/V3xowha9Yk09NDPjToPBBg3gv/3Y3k4LCdRsFNiwbzfZTo1q3k\nbbeZwV4+x3slyflb5rN2v9r82+S//XGh4kh44w3yuutCDhAV9PMPyqiCEdvmK5ZFfvGFadvo0cN0\nx/rI3qN7eceXd7D5x825bs86dxLNzibbtCFfey2kwyzLCvzZB5TfVDBi3XwnPd0M9GrQgJw921fX\nlmXx/fnvMyE5gaNWjXIn0fR0MwZl2rSQDgv6uQeU11QwSoIFwpgxZsj144+bKfQ+snjbYl7w3gX8\nn2//x501SWfMMBHjQ1h5vTSWMhjBu6vjMKIIEfG/9f7uu4HUVKBCBeCii4DPPgOMqHtO82rNsbTr\nUuw+uhtXDr4Sv+75NbIEr7/e9Ah16eL4Gkpzb0lYRKI24RqiQGWj3QJh8WKyeXMz9Do11Te3lmVx\nwMIBjE+O54gVIyJLLDOTbNKE/Pxzx4cE/awDyFtaJSmJFgjZ2SZ0Y3w8+eqrRa5a7jYpO1JY7/16\n7DKxS2TLDy5darqQt293fEjQz9rnfKVVkpJIIMXlMmWAbt3MoKh164CGDYFRo+BHNeXSqpdiSdcl\nOHziMC4fdDnSdqeFl1DTpsCTTwJPPeX4kJycnPB8lTaKUxQANQD8AGA1gFUAnrG3hx1fFVGgsrFm\ngTFrFnnJJWacQ0qKLy4ty+Iniz9hfHI8hy4bGl4ix46ZVcmmTnV8SNDP2Me85F2VBCYU4qX25woA\nfgHQEBHEVw36hsWqBbIcIGmqKR99ZMZuPPmkmcfhA8vTl7PBBw342PjHeDjzcOgJjB9v1s9wOI8m\n6Ofrl9HLKgnJdJLL7M+HAawBcD6A9gCG2rsNBXCH/bkDgJEks0hughGMlsX5UYrHsiyUKVPGf8dl\nypgifloacMYZQKNGQL9+wIkTnrptktgEi7suRraVjZaDWmL1ztWhJdC+PZCUZHp+HODrHJ9YJRR1\ngSkxbAZwNiKIr4ooUNlYt0BZvZps146sU4ccMcIEkvYQy7I4ZOkQxifHc8jSIbRCWYpwzhwzq1VL\nGXnzjndVEp58ySsAWALgDvt72PFVg75hJcUCZ+ZMskULsmlTcsoUz9cUXZWxio3+24gPj32YhzIP\nOT/whhvM5DsHBP1Mfco33vaSiEg5AF8DGEZyvL05Q0Sq2v9PArDT3r4NpqE0l+r2NsVlRARxYaw4\n5RqtW5sFel55xfSs3HgjsHixZ+4uqnIRFnVZhHJx5XDZJ5dhRcYKZwc+9xzw4YeOdjWaoRRKcYoC\nU934HMC7+baHHV8VUaCyJcl8X5inIE6cMA2jSUnkXXd53qPy+bLPGZ8cz48Xf1x8FSU7m6xRg1y2\nzFHaQT9Pr41eVkkAXAPAghGBFNvaIYL4qkHfsJJqUcHhw2Z176pVyQ4dzOhRj1izaw0vHnAxO43p\nxAPHDxS980svka+/7ijdoJ+jD/nE+zYMNy3oG1aSLWo4epR87z0zg/TWW020Ni/cnDjKLhO7sN77\n9bh0+9LCd5wxg2zZ0nG6QT9Hj/OICobaSQtsvEZBHDtG/ve/pkpw443k99970jg6YsUIxifHc8DC\nAQVXUY4dI08/XXtLEJlg6JqeJZggnm2hZGYCI0cCffuaUAHPPQc88IAZ1+ESv+75FfeNuQ/1KtXD\nwNsHouIZFU/doVYtYNYsoE7x64qW5Fms1DADSkFEVaY//XQTL2X5cjPoa8wYoHZt4I03gN9+c8VF\nvcr1MO+JeUgon4BmnzTD4u35emzOOAM4ftwVX6WWSIon4RqioFhWmiyqqih5Wb2a7NaNrFSJvOUW\nM5Q7xHU5C2PUqlFMSE7ge/PfM1UUyyLPOYfcs8fR8UE/My+NWiVRnBDEs3bE0aPA6NHAJ58AGzcC\nnToB998PXHYZEEEpaf3e9eg4piNqVqyJz6p0xTn/090Mb3dAVJXOXIZaJVGcELUvQfnyQOfOwJw5\nwPTpZvXUwGbKAAAHz0lEQVSvBx8ELrwQePVVYOFCIIzp53Ur1cWcx+eg9plJOOf6m7HnLw95cPKl\nCy1hlEJ8D6YUDiSQkgJ89RUwaZIJxNSmDXDTTcA11wB16zorfaxbB3TpguzFiyDp6ShzVgVH7qNW\nXF0gkhKGCkYpJiaEI5ctW4CpU02Us/nzTezY5s2B+vWBCy4ws1LLlzcza/fsATZtAmbPNo2sL74I\nvPCCCeDsEBWMglHBUGJLOHJJTweWLAHWrwc2bDDfjx0DsrOBypWB6tWByy8H2rYFzjwz5ORVMApG\nBUP5nbi4OF2qDqZxONBJfR6jjZ6KK1iW9Xuog0AW6okSSrJYRIreGaVAAlvdK2BKclXEDVQwlEIp\nbaIRteNUoggVDKVISotolPR2C7fQO6QUS0kXDRUL5+hdUhyR2yBa0lCxCA29U0pIBL6OqIuoWISO\n3i0lZEjGfNerikV4FHvHRGSIiGSIyMo82yqJyDQRWSsiU0Xk3Dz/6ykiv4pImoi09erEleDJrabE\n2osXFxcXc+ccLTi5a5/CLPqbl5cBTCNZH8AM+ztEpBGAjjArh7cDMEBE9MmUcHJLHLFAXFycdp9G\ngJNQiT8B2Jdvs4ZJVP5A7ijRaP31VrGInHCfbCLJDPtzBoBE+3M1AFvz7LcVJg6rUorILXFEk3jk\nrhilREbETzPPcmaF7hKpDyV2ySseQTWSagOne4R7FzVMohIyuY2kfgqHNnC6S7h3ciKAzvbnzgDG\n59neSUROE5E6AOoBWBjZKSolDb+EQ9ss3KfYJYhEZCSAVgDiRWQLgDcAvANglIg8AWATgPsAgGSq\niIwCkAogG8BT1CemFELe0aNuL+KjYuENuoCOElXkX8SnTJkyIS/qEytdvEGhK24pJY7clz43fzpZ\nDUwbN52hgqGUGgoTDi1VOEeX6FNKDZZlnVKKiIuLU7HwEefrritKlBBLQ9FLGlrCUBTFMSoYiqI4\nRgVDURTHqGAoiuIYFQxFURyjgqEoimNUMBRFcYwKhqIojlHBUBTFMSoYiqI4RgVDURTHqGAoiuIY\nFQxFURyjgqEoimNUMBRFcYwngiEi7ezYqr+KSA8vfCiK4j+uL9EnImUA/AKgDUxMkkUA7ie5Js8+\nukSfogREtC3R1xLAOpKbSGYB+BIm5qqiKDGOF4JxPoAteb5rfFVFKSF4saan0+rGUg98F0YSgB3q\nT/2pPzSL5GAvBCN/fNUaODWiey4RnXgYJKk/9af+IsOLRs+yMI2eNwDYDhNb9ZRGT0VRYhPXSxgk\ns0WkG4ApAMoAGKxioSglg0AinymKEpvoSE9FURzju2B4MQpURIaISIaIrMyzrZKITBORtSIyVUTO\nzfO/nrb/NBFpG6KvGiLyg4isFpFVIvKMx/7OEJEFIrJMRFJF5N9e+suTRhkRSRGRb7z2JyKbRGSF\n7W+hD/7OFZExIrLGvqeXe/j8GtjXlWsHROQZj6+vp50/V4rICBE53TV/JH0zmDaNdQBqAygHYBmA\nhi6key2ApgBW5tmWDOAl+3MPAO/YnxvZfsvZ57EOQFwIvqoCuNT+XAGmgbehV/7sNMrbf8sCmA/g\nGi/92ek8B2A4gIle3k87jY0AKuXb5qW/oQAez3NPK3p9P+204mC6T2t4mD9rA9gA4HT7+1cAOrvl\nz2/BuBLA93m+vwzgZZfSro1TBSMNQKL9uSqANPtzTwA98uz3PYArIvA7HmYYvOf+AJSHGWp/kZf+\nAFQHMB1AawDfeH0/YQSjcr5tnviDEYcNBWz34/m1BfCTx9dXCeZH7DwYMfwGwI1u+fO7SuLnKNBE\nkhn25wwAifbnajh1XEjY5yAitWFKNgu89CcicSKyzE73B5KrvfQH4F0ALwKw8mzz0h8BTBeRxSLS\nxWN/dQDsEpFPRWSpiAwUkbM89JeXTgBG2p898UdyL4A+AH6DGdawn+Q0t/z5LRiBdMnQSGdRvkM+\nLxGpAOBrAM+SPOSlP5IWyUthfvmvE5HWXvkTkdsA7CSZAqDASUoe3M+rSTYFcDOAp0XkWg/9lYUZ\nNDiAZDMAR2BKul75AwCIyGkAbgcw+g+Jufv86gLoDlPirgaggog85JY/vwXD6ShQN8gQkaoAICJJ\nAHYWcg7V7W2OEZFyMGIxjOR4r/3lQvIAgEkAmnvo7yoA7UVkI8yv4fUiMsxDfyC5w/67C8A4mAmM\nXvnbCmAryUX29zEwApLu8fO7GcAS+xoB767vMgBzSe4hmQ1gLExTgCvX57dgLAZQT0Rq24rbEcBE\nj3xNhGnsgf13fJ7tnUTkNBGpA6AezGhUR4iIABgMIJVkPx/8xee2aIvImTD10RSv/JF8hWQNknVg\nitAzST7s4fWVF5Gz7c9nwdTzV3p4fekAtohIfXtTGwCrYer6rvvLw/04WR3JTdcLf2kArhCRM+28\n2gZAKty6vnAabyIxGKX9BaY1tqdLaY6Eqa+dgGkjeQym8Wc6gLUApgI4N8/+r9j+0wDcFKKva2Dq\n9stgXtwUAO089HcxzES9ZQBWAHiRJxu3XPeXz3crnOwl8er66tjXtgzAqtw84eX1AbgEpvF4Ocwv\ncEWP/Z0FYDeAs/Ns89LfSzAiuBKmR6icW/50pKeiKI7RkZ6KojhGBUNRFMeoYCiK4hgVDEVRHKOC\noSiKY1QwFEVxjAqGoiiO+X9jtWUrEVBaAgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f30b05ffd68>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "# Plot the centroid in red\n",
    "ax.scatter(C_x, C_y, color='red', marker='+', s=100)\n",
    "\n",
    "# Plot the major axis\n",
    "x1 = C_x + semi_major_length * np.cos(theta)\n",
    "y1 = C_y + semi_major_length * np.sin(theta)\n",
    "x2 = C_x - semi_major_length * np.cos(theta)\n",
    "y2 = C_y - semi_major_length * np.sin(theta)\n",
    "ax.plot([x1, x2], [y1, y2], color='green', lw=1)\n",
    "\n",
    "# Plot the minor axis\n",
    "x1 = C_x + semi_minor_length * np.cos(theta + np.pi/2)\n",
    "y1 = C_y + semi_minor_length * np.sin(theta + np.pi/2)\n",
    "x2 = C_x - semi_minor_length * np.cos(theta + np.pi/2)\n",
    "y2 = C_y - semi_minor_length * np.sin(theta + np.pi/2)\n",
    "ax.plot([x1, x2], [y1, y2], color='green', lw=1)\n",
    "\n",
    "# Plot the ellipse\n",
    "angles = np.arange(0, 360, 1) * np.pi / 180\n",
    "x = 0.5 * major_length * np.cos(angles)\n",
    "y = 0.5 * minor_length * np.sin(angles)\n",
    "\n",
    "R = np.array([[np.cos(theta), -np.sin(theta)],\n",
    "             [np.sin(theta), np.cos(theta)]])\n",
    "\n",
    "x, y = np.dot(R, np.array([x, y]))\n",
    "x += C_x\n",
    "y += C_y\n",
    "ax.plot(x, y, lw=1, color='red')\n",
    "\n",
    "# Show the image\n",
    "_ = ax.imshow(a, interpolation=\"none\", cmap=\"gray\", origin='lower', aspect='equal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems to work but there is an underestimation of the axes length...\n",
    "\n",
    "The method is not exactly the same but this understimation does not appear in the [ImageJ/Fiji algoritm](http://imagej.nih.gov/ij/source/ij/process/EllipseFitter.java).\n",
    "\n",
    "Any modifications are welcome !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "Some links which help me : \n",
    "\n",
    "- Kieran F. Mulchrone, Kingshuk Roy Choudhury, Fitting an ellipse to an arbitrary shape: implications for strain analysis, Journal of Structural Geology, Volume 26, Issue 1, January 2004, Pages 143-153, ISSN 0191-8141. http://dx.doi.org/10.1016/S0191-8141(03)00093-2\n",
    "- https://en.wikipedia.org/wiki/Image_moment\n",
    "- https://github.com/shackenberg/Image-Moments-in-Python/blob/master/moments.py\n",
    "- http://imagej.nih.gov/ij/source/ij/process/EllipseFitter.java\n",
    "- http://imagej.1557.x6.nabble.com/Ellipse-fitting-algorithm-used-by-ImageJ-td3686611.html\n",
    "\n",
    "Another method to detect ellipse : \n",
    "\n",
    "- Xie, Y. X. Y., & Ji, Q. J. Q. (2002). A new efficient ellipse detection method. Object Recognition Supported by User Interaction for Service Robots, 2(c), 0–3. http://doi.org/10.1109/ICPR.2002.1048464"
   ]
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}