RANSAC_en.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RANSAC\n",
    "\n",
    "It is well known that least squares parameter estimation is not resistent, i.e. sensitive to outliers. Even one outlier may break least squares estimation, which means a possibly infinitely large parameter bias.\n",
    "\n",
    "This problem is effectively solved by **RANSAC (RANdom Sample And Consensus)** introduced By Fischler and Bolles in 1981. They solved a problem with a large percent of outliers in the input data. Even 50% of outliers may be tolerated by RANSAC.\n",
    "\n",
    "Therefore RANSAC is a resistent iterative parameter estimation procedure that uses a subset of data consistent with a model. At each iteration step the procedure works as follows:\n",
    "\n",
    "1. Select a random minimal sample from data and check conformity\n",
    "2. Estimate model and check its validity\n",
    "3. Classify data as conformal or outlier based on deviations from the model - data are conformal if deviations are below threshold\n",
    "4. Keep model if the number of conformal data is maximum.\n",
    "\n",
    "Iteration ends upon reaching a specified maximum number or upon a specific criterion. Model parameters are then estimated using the best model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RANSAC ellipse fitting\n",
    "\n",
    "RANSAC is illustrated with a simple problem.\n",
    "\n",
    "Generate conformal data that approximately fit on an ellipse:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "t = np.linspace(0, 2 * np.pi, 50)\n",
    "a = 10\n",
    "b = 5\n",
    "xc = 20\n",
    "yc = 30\n",
    "theta = np.pi/6\n",
    "x = xc + a*np.cos(theta)*np.cos(t) - b*np.sin(theta)*np.sin(t)\n",
    "y = yc + a*np.sin(theta)*np.cos(t) + b*np.cos(theta)*np.sin(t)\n",
    "data = np.column_stack([x, y])\n",
    "# for reproducibility:\n",
    "np.random.seed(seed=1234)\n",
    "data += np.random.normal(size=data.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add some outliers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAUGUlEQVR4nO3de4yldX3H8feXGRYKDeE2yGWhi5VoFbSSiWW8tKOLESwBGrDBmLoB4lql9UINuCUWG9LgLRWJlboRZGkoF1mUS6pItk5NkxEyKwjIRbaCsAIyEqCNFxZ2v/3jeY5zOHtmZ2fOmTnn/Pb9SiZnnss557sPM5/58X1+53kiM5EklWW3XhcgSeo+w12SCmS4S1KBDHdJKpDhLkkFGu51AQAHHnhgrlixotdlSNJA2bhx4y8zc6Tdtr4I9xUrVjA1NdXrMiRpoETEz2bbZltGkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwlzTYJifh4ourx0GziLX3xTx3SVqQyUlYuRK2bIFly2DDBhgb63VVO2eRa3fkLmlwTUxU4bh1a/U4MdHrinbeItduuEsaXOPj1ah3aKh6HB/vdUU7b5Frty0jaXCNjVXtjImJKhwHpSUDi1579MNt9kZHR9Nry0jS/ETExswcbbfNtowkFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklSgOcM9Iq6IiKcj4r6mdZ+PiAcj4p6I+GZE7Nu0bU1EbIqIhyLiXYtVuCRpdjszcr8SOKFl3e3A0Zn5euAnwBqAiHgtcAbwuvo5X4mIoa5VK2nwDPKdkgbYnJf8zczvR8SKlnXfbVr8AXB6/f0pwLWZ+QLwSERsAt4E+F9V2hUN8p2SBlw3eu5nAd+uvz8MeLxp2+Z63XYiYnVETEXE1PT0dBfKkNR3BvlOSQOuo3CPiAuAl4CrG6va7Nb2gvGZuTYzRzNzdGRkpJMyJPWrQb5T0oBb8J2YImIVcBKwMmfu+LEZOLxpt+XAEwsvT9JAG+Q7JQ24BYV7RJwAnA/8WWb+umnTzcC/R8Q/A4cCRwF3dlylpME1Nmao98Cc4R4R1wDjwIERsRm4kGp2zB7A7REB8IPM/OvM/HFEXA/cT9WuOSczty5W8ZKk9ryHqiQNKO+hKkm7GMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSrQnOEeEVdExNMRcV/Tuv0j4vaIeLh+3K9eHxFxaURsioh7IuLYxSxektTezozcrwROaFn3SWBDZh4FbKiXAU4Ejqq/VgOXdadMSYtqchIuvrh6VBGG59ohM78fEStaVp8CjNffrwMmgPPr9VdlZgI/iIh9I+KQzHyyWwVL6rLJSVi5ErZsgWXLYMMGGBvrdVXq0EJ77q9oBHb9eFC9/jDg8ab9NtfrthMRqyNiKiKmpqenF1iGpI5NTFTBvnVr9Tgx0euK1AXdPqEabdZlux0zc21mjmbm6MjISJfLkLTTxserEfvQUPU4Pt7ritQFc7ZlZvGLRrslIg4Bnq7XbwYOb9pvOfBEJwVKWmRjY1UrZmKiCnZbMkVYaLjfDKwCPlM/3tS0/m8i4lrgT4Dn7bdLA2BszFAvzJzhHhHXUJ08PTAiNgMXUoX69RFxNvAY8J569/8A3g1sAn4NnLkINUuS5rAzs2XeO8umlW32TeCcTouSJHXGT6hKUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQXqKNwj4uMR8eOIuC8iromIPSPiyIi4IyIejojrImJZt4qVJO2cBYd7RBwGfAQYzcyjgSHgDOCzwBcz8yjgWeDsbhQqSdp5nbZlhoHfi4hhYC/gSeAdwA319nXAqR2+hyRpnhYc7pn5c+ALwGNUof48sBF4LjNfqnfbDBzW7vkRsToipiJianp6eqFlSJLa6KQtsx9wCnAkcCiwN3Bim12z3fMzc21mjmbm6MjIyELLkCS10Ulb5njgkcyczswXgRuBNwP71m0agOXAEx3WKEmap07C/THguIjYKyICWAncD3wPOL3eZxVwU2clSpLmq5Oe+x1UJ05/CNxbv9Za4Hzg3IjYBBwAXN6FOiVJ8zA89y6zy8wLgQtbVv8UeFMnrytJ6oyfUJWkAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCtRRuEfEvhFxQ0Q8GBEPRMRYROwfEbdHxMP1437dKlaStHM6Hbl/CfhOZr4GeAPwAPBJYENmHgVsqJclSUtoweEeEfsAfwpcDpCZWzLzOeAUYF292zrg1E6LlCTNTycj91cC08DXI+KuiPhaROwNvCIznwSoHw9q9+SIWB0RUxExNT093UEZkqRWnYT7MHAscFlmvhH4FfNowWTm2swczczRkZGRDsqQJLXqJNw3A5sz8456+QaqsP9FRBwCUD8+3VmJkqT5WnC4Z+ZTwOMR8ep61UrgfuBmYFW9bhVwU0cVSpLmbbjD5/8tcHVELAN+CpxJ9Qfj+og4G3gMeE+H7yFJmqeOwj0z7wZG22xa2cnrSpI64ydUJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVKCOwz0ihiLiroi4tV4+MiLuiIiHI+K6iFjWeZmSpPnoxsj9o8ADTcufBb6YmUcBzwJnd+E9JEnz0FG4R8Ry4M+Br9XLAbwDuKHeZR1waifvIUmav05H7pcA5wHb6uUDgOcy86V6eTNwWLsnRsTqiJiKiKnp6ekOy5AkNVtwuEfEScDTmbmxeXWbXbPd8zNzbWaOZuboyMjIQsuQJLUx3MFz3wKcHBHvBvYE9qEaye8bEcP16H058ETnZQ6YyUmYmIDxcRgb63U1knZBCw73zFwDrAGIiHHgE5n5voj4BnA6cC2wCripC3UOjslJWLkStmyBZctgwwYDXtKSW4x57ucD50bEJqoe/OWL8B69MzkJH/pQ9TU5uf32iYkq2LdurR4nJpa6QknqqC3zO5k5AUzU3/8UeFM3XrfvTE7C298OL7xQLV9xRRXezSPz8fFqxN4YuY+P96BQSbs6P6HabHISLr64/YgcZkblDS++uP3IfGysasVcdJEtGUk905WRexGae+XDw3DmmfD+97cflTdG7rvv3n5kPja2fah7klXSEtq1wn1HAdvcK9+6Fb76VVi37uWj77ExuPRSuPxyOPRQOO+8mW2N1z7gAHjmmZnQb6z72Mc8ySppyew64T7XLJbxcRgaqoIdIHPmhGhzgDdC+t57q3Bvfu0XXoBt22C33arRfwS89FK1vHVrtW3LFrjqKkfxkhbVrhPu7WaxtAZrRPWVWQVy6wnR2V6jsX5b/UHdbduqfjxUr9V4vYgq9K+4onoNR/GSFkn5J1QbUxfvvLMamQ8NtZ/FMjFRjbIbQXz88e1H98uWbf8ajfW71Ydzt92qfnxj3z32gHPPrUb3J5440/pxqqSkRVL2yH1ysgrexgyX3XeHD3xg5kRpcw++dQrjpz+9/Yi6MROmtaXSvL61537VVfDUU3DJJVWgDw1Vo3dwqqSkRVNeuDcH9sTETHsEqpH5EUfMBHtrD75dcLdqNxNmtvWTk9VJ2d/+tvo/AqhaNh/8YFWHPXdJi6SscG8N7EsuqUbrjZF780i5Xf98zZruhm3jPbLp2mmZsM8+1XtJ0iIpK9xbA/uuu+Css6q2yMEHv3ze+lJ8krTxHr/5zcvX3313999LkpqUFe7NgT3XrJTZ+ufd1HiPz30OvvWtmfWnndb995KkJmWF+9hY1YpZvx722gtuuWXHUx9n6593S6P/f9551SyZ9eurYF+9evHeU5IoLdybP2Q0PFzNTIHezEpp/mDT0BB8+ctw221LW4OkXVZZ4d7cc4dq2mOvZqVMTMx8YnXbNjjnHDjmGGfHSFoSZYR783Vdmk+Stl74qxvvsbN/KBqXM2j+1Gq71pAkLYLBD/d20x8bHyLqZrDPdXel1vAfG6taMeecUwX7Hnv4gSVJS2bww711+uMzz3R/Dvlc16WZLfxXr65aMV4kTNISG/xwP+CA6oJc7S701S1zzYnfUfgv9owcSWpjsMO9MTumcQL1bW+buRBXNwO1MV/9qqu2f/92vX7bL5J6bLDDfWLi5ddt+e534fbbYc89X94Xn+/J0LVr289JX7euCvB166refvMNOBaj1y9JCzTY4T4+PnP99YbWm2zszMnQZmvXVhf2guqPBVQB39p6Wb++fa+/cR9WQ15SDw3+9dzf+taXL0dUUxB3dIGwHVm/vv1y67XcTztt+2u7N/6QfOpT1eNsN9qWpEU2uCP35hH57rvDq14FDz88cyOMe++t9rvzzuqx+YTrjto0p502M2JvLEP7a9G0zoS5+OK57/YkSUtgcMO99dOohx8ODz1UtWW2boUPf7gK9Mb13IeGqr447LhN0+ixN26CfcwxM9taZ760Li/FlSYlaScMblumXZukcS0ZePl9TBvLzzyzc22aY46pRv633DK/9kpjdH/RRd4bVVJPDe7IfbZL9jY+ETo8XI3iGwHfPJKea3S9MzfT3lFdhrqkHhvccIftg7T1E6EwMze9+TozrX8UWnvwtlckDbjI5mmE83lixOHAVcDBwDZgbWZ+KSL2B64DVgCPAn+Zmc/u6LVGR0dzampqQXV0bLapkvOdGy9JSywiNmbmaLttnfTcXwL+LjP/CDgOOCciXgt8EtiQmUcBG+rl/jVbD35srPv3VJWkJbLgcM/MJzPzh/X3/wc8ABwGnAKsq3dbB5zaaZGLqvXErC0YSQXoSs89IlYAbwTuAF6RmU9C9QcgIg6a5TmrgdUARxxxRDfKWJiluJeqJC2xBffcf/cCEb8P/BfwT5l5Y0Q8l5n7Nm1/NjP329Fr9LTnLkkDarF67kTE7sB64OrMvLFe/YuIOKTefgjwdCfvIUmavwWHe0QEcDnwQGb+c9Omm4FV9fergJsWXp4kaSE66bm/Bfgr4N6IuLte9/fAZ4DrI+Js4DHgPZ2VKEmarwWHe2b+NxCzbF650NeVJHVucK8tI0maleEuSQXqeCpkV4qImAZ+1us62jgQ+GWvi1gga+8Na++NXbX2P8jMkXYb+iLc+1VETM02h7TfWXtvWHtvWPv2bMtIUoEMd0kqkOG+Y2t7XUAHrL03rL03rL2FPXdJKpAjd0kqkOEuSQUy3KluGRgR34uIByLixxHx0Xr9/hFxe0Q8XD/u8NLFvRQRQxFxV0TcWi8fGRF31LVfFxHLel1jOxGxb0TcEBEP1sd/bFCOe0R8vP55uS8iromIPfv1uEfEFRHxdETc17Su7XGOyqURsSki7omIY3tX+ay1f77+mbknIr4ZEc2XGV9T1/5QRLyrN1X/rpbtam/a9omIyIg4sF7u6nE33Csl3DLwo1R3w2r4LPDFuvZngbN7UtXcvgR8JzNfA7yB6t/Q98c9Ig4DPgKMZubRwBBwBv173K8ETmhZN9txPhE4qv5aDVy2RDXO5kq2r/124OjMfD3wE2ANQP17ewbwuvo5X4mIoaUrdTtXsn3tjXtQv5Pq4ooN3T3umelXyxfVZYrfCTwEHFKvOwR4qNe1zVLvcqpfzncAt1Jd0O2XwHC9fQy4rdd1tql7H+AR6hP7Tev7/rhT3VLycWB/qgvw3Qq8q5+PO9VN6++b6zgDXwXe226/fqm9ZdtfUN1TAqqQX9O07TZgrN9qB26gGsw8Chy4GMfdkXuLHd0yEGh7y8A+cAlwHrCtXj4AeC4zX6qXN1OFUb95JTANfL1uKX0tIvZmAI57Zv4c+ALVyOtJ4HlgI4Nx3BtmO86NP1wN/f7vOAv4dv1939ceEScDP8/MH7Vs6mrthnuT+paB64GPZeb/9rqenRERJwFPZ+bG5tVtdu3HOa/DwLHAZZn5RuBX9GELpp26P30KcCRwKLA31f9Wt+rH4z6XQfn5ISIuoGqrXt1Y1Wa3vqk9IvYCLgD+od3mNusWXLvhXhvgWwa+BTg5Ih4FrqVqzVwC7BsRjev1Lwee6E15O7QZ2JyZd9TLN1CF/SAc9+OBRzJzOjNfBG4E3sxgHPeG2Y7zZuDwpv368t8REauAk4D3Zd3HoP9r/0OqAcGP6t/Z5cAPI+Jguly74c5g3zIwM9dk5vLMXEF1Iuk/M/N9wPeA0+vd+rX2p4DHI+LV9aqVwP0MwHGnasccFxF71T8/jdr7/rg3me043wy8v569cRzwfKN90y8i4gTgfODkzPx106abgTMiYo+IOJLq5OSdvaixncy8NzMPyswV9e/sZuDY+nehu8e9lyca+uULeCvV//7cA9xdf72bqne9AXi4fty/17XO8e8YB26tv38l1Q/1JuAbwB69rm+Wmv8YmKqP/beA/QbluAP/CDwI3Af8G7BHvx534BqqcwMv1oFy9mzHmao98C/A/wD3Us0I6rfaN1H1pxu/r//atP8Fde0PASf2W+0t2x9l5oRqV4+7lx+QpALZlpGkAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUD/DxTlOXCPOQvVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data[0] = (100, 100)\n",
    "data[1] = (110, 120)\n",
    "data[2] = (120, 130)\n",
    "data[3] = (140, 130)\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(data[:,0],data[:,1],'r.')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use total least squares estimation\n",
    "\n",
    "The image processing library [Scikit-image](http://scikit-image.org) will be used for estimation. It contains implementation of both least squares and RANSAC estimation that can be used for our problem. We import only a few functions in `sm.py`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[73.53514375565935,\n",
       " 78.26758207018908,\n",
       " 9.294963688237154,\n",
       " 80.77273222385621,\n",
       " 2.3039489382317697]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import sm\n",
    "model = sm.EllipseModel()\n",
    "model.estimate(data)\n",
    "# ellipse parameters:\n",
    "# xc, yc, a, b, theta\n",
    "np.set_printoptions(suppress=True)\n",
    "model.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we see these parameters are completely wrong. Original semi-axes of the ellipse were 10 and 5, the centre was at $O(20,30)$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RANSAC estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[20.060780840092757, 29.44604087550921, 5.615422139522865, 9.2913820935979, 2.092276704111832]\n",
      "[-0.06078084  0.55395912  4.38457786 -4.29138209 -1.56867793]\n"
     ]
    }
   ],
   "source": [
    "# at least 5 points are needed for ellipse fitting\n",
    "n_min = 5\n",
    "# maximum distance of conform points\n",
    "t_max = 3.0\n",
    "ransac_model, inliers = sm.ransac(data, sm.EllipseModel, n_min, t_max, max_trials=50)\n",
    "print(ransac_model.params)\n",
    "original_params = np.array([xc,yc,a,b,theta])\n",
    "print(original_params - ransac_model.params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Conformal data were found and stored in array 'inliers' as logical 'True' elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([False, False, False, False,  True,  True,  True,  True,  True,\n",
       "       False,  True,  True,  True,  True,  True,  True,  True,  True,\n",
       "        True,  True,  True,  True,  True,  True,  True,  True,  True,\n",
       "        True,  True,  True,  True,  True,  True,  True,  True,  True,\n",
       "        True,  True,  True,  True,  True,  True,  True,  True,  True,\n",
       "        True,  True,  True,  True,  True])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inliers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the first four data are outliers that were added by us."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot best fitting ellipse to conformal data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXhU5fXA8e8hsgqKrFYWwRUtVdAIIrYKLqDYqlXrgtZqlWLda/0htVbFqmjVlrqjoFipiAug4gZiQCSiQWVfFdcgi4qyBpKc3x/npsQ4SSbJzNyZO+fzPPNMcmducoYhZ977LucVVcU551x01Qs7AOecc8nlid455yLOE71zzkWcJ3rnnIs4T/TOORdxO4UdQCytWrXSTp06hR2Gc85ljDlz5qxT1daxHkvLRN+pUycKCgrCDsM55zKGiHxa2WPedeOccxHnid455yLOE71zzkWcJ3rnnIs4T/TOORdxnuidcy7iqp1eKSKNgBlAw+D5z6rqjSLyFtAseFob4F1VPSXG+SXA/ODbz1T1VwmJ3DmX/vLzIS8Pjj4aevUKO5qsFc88+iKgr6puFJH6wEwReUVVf172BBF5DphUyflbVLVbAmJ1zmWS/Hw45hjYtg0aNIA33vBkH5Jqu27UbAy+rR/c/lfEXkSaAX2BiUmJ0DmXmfLyLMmXlNh9Xl7YEWWtuProRSRHRD4E1gBTVHV2uYdPBd5Q1e8rOb2RiBSIyDsi8qOunXK/Y1DwvIK1a9fG/QKcc2nq6KOtJZ+TY/dHHx12RFkrrhIIqloCdBOR5sAEEemqqguCh88GHq3i9I6qWigiewHTRGS+qn4U43eMBEYC5Obm+rZXzmW6Xr2su8b76ENXo1o3qrpeRPKA/sACEWkJ9MBa9ZWdUxjcfxyc2x34UaJ3zkVQr16e4NNAtV03ItI6aMkjIo2BY4ElwcNnAC+p6tZKzt1NRBoGX7cCegOLEhG4c865+MTTR/8T4E0RmQe8h/XRvxQ8dhbwVPkni0iuiJR15RwAFIjIXOBNYLiqeqJ3zrkUEtX06w7Pzc1VL1PsnHPxE5E5qpob6zFfGeuccxHnid455yLOE71zzkWcJ3rnMkF+Ptx+u907V0NpuWesc66cWDVjwBciubh5oncu3VWsGfPEEzBmjBcLc3Hzrhvn0l3FmjHgxcJcjXiL3rl0V7FmDPywRe/Fwlw1PNE7lwkq1ozxYmGuBjzRO5dBSkpg/XpY16IXm/r1oqgIivKw+yJr5NerB40aQcOGdl92a9ECWrWC+vXDfhUu1TzRO5cmNm6ETz6BlSt/eL9mDaxbZ7dvvoG6Vi1p2RLatLHbT34C++4L++0H++9v97vumoAX49KKJ3rnUmzDBli4EObNg/nz7bZoEVTcb6dxY+jUyZLxwQdba7zs1rIlNG1qrfYGDXbcN2hgHwRbt1oLv+x+yxb7kFi92j441qyxr999F8aPh9LSHb+3bVvo1g0OP9x6hXr0gN12S+k/kUswT/TOJdHWrfD++zYVPj/fvl65csfjTZtC167wq1/B3ntD5847bq1bg0jyYywqgo8/hqVLYdkyWLIECgpg2LAdVw9dusBRR8GAAdC3L+y8c/Ljconj1SudS6C1a22MdNasHYl9+3Z7bK+9IDcXDjoIfvYzu+25p/Wpp6MNG+C99+Cdd+y15OVZ91LDhtCnD5x0EpxyCrRrF3akDqquXumJ3rk62LoVZs6EKVPs9sEHdrxxY0vqvXrBEUdYN0jbtuHGWlfbtsFbb8FLL8HkybB8uX1IHXccXHABnHyyDfq6cHiid66isiZqLaYnrlwJEyfCK69Y4tu61Way9OplSe/YY+HQQ6M/u2XpUhg7Fh5/HD7/3Prxzz4brrjCBnZdanmid668WLVjqkj2qjZwOnEiTJgAc+fa8QMPtMR+3HHWf920aYriTzMlJTBtGjz2GDz/vP2znnYaDB0KhxwSdnTZo6pE74OxLvtUrB2Tl/ejRK9q/ev//a8l95UrbWC0d2+46y7rm95771CiTzs5OTs+8NasgX/9C+6/H559Fvr1g5tusq4rF540HQZyLokq1o4pV0Lgk0/g1luttZ6bC/fdBwccACNHwqpV1lVzzTWe5CvTpg3cdht89pndf/CBfYaedx4UFoYdXfbyrhuXncr10a8/oBfjx8OTT1oiB/jFL+Dcc+H0030OeV1s3Ghl9O+6y8YsbrgBrrrKZu64xPI+eucqULXFQg89BOPG2YBqly7W8jznHFuo5BLno4/sSmjSJLtaeuopm2bqEqdOm4OLSCMReVdE5orIQhG5OTj+uIisFJEPg1u3Ss4/X0SWB7fz6/ZSnKubDRvg4YdtkPDww+GZZ+D8822++KJF8Je/eJJPhr33tsHsyZNthe5hh1lffvkVuS554hmMLQL6qupGEakPzBSRV4LHrlXVZys7UURaADcCuYACc0TkBVX9tq6BO1cTS5fCiBHWPbNhg7UmH3zQWu+77BJ2dNnjxBNtBtNFF8HVV8Orr9o+Km3ahB1ZtFXbolezMfi2fnCLt7+nHzBFVb8JkvsUoH+tInWuhlRhxgwrL9ClC4weDb/+tXXPf/ghDB5cTZL3fVqTonVra90/+CBMn25XVkuW1PKH+XsUl7hm3YhIjoh8CKzBEvfs4KFbRWSeiPxTRGINr7QDPi/3/RfBsVi/Y5CIFIhIwdqK1Z2cq4HiYivU1bOnzW/Pz4cbb7SZII8/boml2hoyZXPtb7jB7j2RJJSIfdBOnw6bNtnMnOnTa/hD/D2KW1yJXlVLVLUb0B7oISJdgaFAF+AwoAUwJMapsf6cYl4NqOpIVc1V1dzWrVvHFbxz5RUVwQMPWNndM8+0uu0PPgiffmpzuWvUPRBrrr1LuB49rJZO27Zw/PE2MB43f4/iVqN59Kq6HsgD+qvqqqBbpwh4DOgR45QvgA7lvm8P+Gxal1BlCX7vveHSS62s74QJsHixtRqbNKnFD61irr1LrM6drQhcz54wcCA891ycJ/p7FLd4Zt20FpHmwdeNgWOBJSLyk+CYAKcAC2Kc/hpwvIjsJiK7AccHx5yrs4oJvlMnmDoV3n7bVq7m5NThh5ft03rLLdWWSMgKSe4Lb9HCagf17Gn1cl6LJ0v4exQ/Va3yBhwEfADMw5L534Lj04D5wbEngabB8Vzg0XLnXwisCG4XVPf7VJVDDz1UnavMtm2qDzyg2q6dKqj27q06dapqaWnYkUXUrFmqjRur5uTY/axZSftV336revDB9mtmzkzar4kkoEAryanVTq9U1XlA9xjH+1by/ALgonLfjwZGx/ex41zlVK1L5rrrrERu794wZoxthJGKDTqyVhy1gRKl+eJ8Xh/wLj//dhAnn9yYDz6ADh2qP89VzWvduIwwaxYceaRVRaxf32qiv/WWTbbwJJ9kqeoLD2bRtLnjGl5c3ZNtW4r5zW/ss8XVjSd6l9amTrXk3ru3VZB85BErEzxggCf4lElVX3i5K4f9ihcx6uQXeecdGBJrPp+rES9T7BKjDht5xLJhg5W4LRv7u/lmq5Xie5WGpFev+N/X2v5fKLtyCPYJOOPy3bmitZVK6NPHFr652vGiZq7uariRR1VUbbHTNdfAl1/asfnzbQNtlwHq+n+hwofEtm1WLnr9epsu6x/0latTUTPnqpWghSuLF9s2fGedZYubZs2yxJ+UJO9L55Ojrv8XevWyramCD4cGDWzR2+efw7BhCY82a3iid3VXx8G6TZusH/agg2xXpwcesGqSSZsW7UvnkycJA7e9e8OFF8I998DChXX+cVnJE72ruzoM1k2bZgn+zjvht7+FZcvgkkvquNipOr50PnmSNHB7xx3QrBlce21CflzW8cFYlxg1GawDvvvO/mgfeQT22ccKWv3iF0mMr7wKg36+dD7Bavh/IR6tWsGf/wzXX2+VR7vF3P3CVcZb9C7lJk+Gn/4URo2yZD9vXgqTPPjS+Qz1xz9aq3748LAjyTye6F3KfPutbdV30knQvLl1jd95JzRuHEIwFQb9XPpr3tyS/TPP2MroOsuiAXlP9C4l8vKsL37cOPjb32DOHCtR61xNXH017LQT3H9/HX9Qlg3Ie6J3SbVtmzWc+/a1lnt+vi1+ahhrmxrnqtG2rV0RjhtnY+m1lmUD8p7oXdIsXQpHHGF9qr//vU2dzI25nMO5+J1zDqxebTO2ai3Latl7oncJpwqPPgqHHGL1aZ57zmbXNG0admQuCgYMsL1+//vfOvyQLBuQ9+mVLqE2b7Zdnf7zH+v6HDMG2sXcJdi52mnUyDZ5nzDBGhS1XnORhGmg6cpb9C5hli+3jbeffNL64V9/3ZO8S47jjrO1GHPnhh1JZvBE7xLi+eet/72w0LaE+9vfoJ7/73JJ8vOf2/2MGeHGkSn8T9HVyfbttmLxtNOgSxcbcO3XL+yoXNR16GCbinuij48neldr33wD/fvD3Xfb5twzZkDHjmFH5bLFL35hu4ylYaX1tOOJ3tXKsmXWHz9zJjz+ONx3n8+Nd6nVrRusWwffvPZe1qxwrS2fdeN2iHNnoDfegNNPt71bp02zMrLOpdo++9j9ilP+TMvit+u86U2UVduiF5FGIvKuiMwVkYUicnNwfKyILBWRBSIyWkTqV3J+iYh8GNxeSPQLcAkS55Lwhx+2Pvj27eHddz3Ju/D8L9Fv65g1K1xrK56umyKgr6oeDHQD+ovI4cBYoAvwM6AxcFEl529R1W7BzXd9TFfVLAkvfTufq498l8GDrV/+7behU6cwAnXOdO4MIsqKnP2zZoVrbVXbdaO2qezG4Nv6wU1V9eWy54jIu0D7pEToUqOKGu1F09/h/GO+4OmSM7gq517uui6XnF388tilQBXdiQ0bwu67C5/nXgy9chK2MX0UxdVHLyI5wBxgH+B+VZ1d7rH6wHnAlZWc3khECoBiYLiqTqxbyC4pypaEV/ij2rABfn1RW6aWHM6dXMu1/BPeugWO9D8ol2RxbDTerBlsbNLWKue5SsWV6FW1BOgmIs2BCSLSVVUXBA8/AMxQ1bcqOb2jqhaKyF7ANBGZr6ofVXySiAwCBgF09Dl64aiwJHzNGqsr8sHKTjxe/2LOL33ML49d6sTqTqyQ6Js2hY0bY57tyqnRrBtVXS8ieUB/YIGI3Ai0Bv5QxTmFwf3HwbndgR8lelUdCYwEyM3N9ZmxIVu50gZdv/gCJk0SBrS4EPL28stjlzpxbPnoiT4+1SZ6EWkNbA+SfGPgWOAOEbkI6Acco6qllZy7G7BZVYtEpBXQG7gzceG7ZFi2DPr0gS1bYOpUKzUM2VMAyqWJSroTy2va1EoWu6rF06L/CTAm6KevB4xX1ZdEpBj4FMgXEYDnVXWYiOQCg1X1IuAA4GERKQ3OHa6qi5LySlxCLFlim4QUF9tK165dw46IuOf3uwiqpsJkSUkdqldmkXhm3czDulsqHo95rqoWEEy1VNVZ2PRLlwEWL7Ykr2p59cADw46IuAbkXPbatAl23jnsKNKfl0BwACxaZN01AG++mSZJHrJuyzdXM5s2QZMmYUeR/jzROxYutF6RevUsjx5wQNgRlZNlW765mvlBiz4/32veVMJr3WS5jz6ynpH69a0lv99+YUdUQRwDci57bdgQbFHpXXxV8kSfxQoLbaee4uI0TfJlsmjLNxe/bdvgq6+s7lI8c+6zmSf6LPXNNzZPfu1aq0CZVt01zsXh889t4kDnzsC2ltb3qOpdfDF4os9CmzbBSSfZfPmXX4bDDgs7Iudq7pNP7L7TpoVw7VXWmq9XD/71L2/NV+CDsZmuhgNQ27bZtn+zZ8O4cdat6VwmWrnS7jt9Ot3+Y5eWWov+66/DDSwNeYs+k9VwAEoV/vhHeO01GDUKTj01hbE6l2DLl9skgva/7A73VV0qIdt5os9kNRyAuusuS/B//StceGHKonQuKebMgYMOgp1+7jOzquOJPpPFUfSpzMSJMGQInHEG3HxzyiJ0LilKS6GgAM4+OzjgM7Oq5Ik+k8U5x/z992HgQBt0HTPGxqucy2QrVsB330FubtiRZAZP9JmumpZMYSH88pfQsiVMmgSNG6cwNueS5L337N5njMXHE32Ebd8OZ54J69fDrFmw++5hR+RcYrz9tpU+SJuaTGnOE32EDR0KM2fC2LFw8MFhR+NcYqjCq69apdWdPIPFxXtrI+r55+Huu2065TnnhB2Nc4mzbJnNoT/hhLAjyRye6CNo+XK44ALo0QPuuSfsaJxLrFdesXtP9PHzRB8xW7bA6afbJe348dCwYRVP9rKuLgO98gp06QKdOoUdSebwHq6IGToU5s2zGjZ77lnFE72sq8tA69fD9OnWJeni5y36CHnjDRgxAi6/PI7LWt+5yWWg556DoqJyC6VcXDzRR8T69fC738H++8Pw4XGc4Ds3uQz05JO2b4IvlKoZ77qJiMsug1WrrEcmrj00fecml2E++8z+uw4bBiJhR5NZPNFHwDPP2Fz5m26q4UpBrw/iMshTT9n9wIHhxpGJqu26EZFGIvKuiMwVkYUicnNwvLOIzBaR5SLytIg0qOT8oSKyQkSWiki/RL+AbPfNN3DppXYp+5e/hB2Nc8lRWgqPPQZHHAF77RV2NJknnj76IqCvqh4MdAP6i8jhwB3AP1V1X+Bb4PcVTxSRA4GzgJ8C/YEHRCQnUcE7m2Xz9dfwyCNWm9u5KHr9dVi6FC65JOxIMlO1iV7NxuDb+sFNgb7As8HxMcApMU4/GRinqkWquhJYAfSoc9QOsPo1I0fClVdCt25hR+Nc8owYYbWafvObsCPJTHHNuhGRHBH5EFgDTAE+AtaranHwlC+AdjFObQd8Xu77yp6HiAwSkQIRKVi7dm288Wet7dth8GBo397ry7toW7LEattccolNEHM1F1eiV9USVe0GtMda5AfEelqMY7HGxmM9D1Udqaq5qprbunXreMLKaiNGwPz58O9/Q7NmYUfjXPLce68l+D/8IexIMleN5tGr6nogDzgcaC4iZbN22gOFMU75AuhQ7vvKnudqYM0aa8UPGACnxOowcy4i1q2Dxx+3BVJt24YdTeaKZ9ZNaxFpHnzdGDgWWAy8CZwePO18YFKM018AzhKRhiLSGdgXeDcRgWezYcNgyxbl7v0eRt7xOjXJNHas1VSpV8/ux44NO6Ls8o9/WP2mIUPCjiTDqWqVN+Ag4ANgHrAA+FtwfC8saa8AngEaBsd/BQwrd/71WJ/+UuCE6n6fqnLooYeqi23ZMtWdckp0cM5I1Zwc1caNVWfNCjusSHrySdUmTVStArrdmjSx4y75vvrK/r0HDgw7kswAFGglOVXs8fSSm5urBQUFYYeRls44A155YRsrijuze2mhlTC45RabZ+kSqlMn+PTTHx/fc0/45JNUR5N9/vQnG4tavNjKHriqicgcVY1ZHMJr3WSQ2bPh2Wfhz+euZveG33qdmiT77LOaHXd1VK5sdmEhPPggnHeeJ/lE8BIIGeT666FNG7jmXx3gIq9Tk2wdO8Zu0XfsmPpYIq9C2exbT/iY7dt354Ybwg4sGjzRZ4j33rMaZHfeGUyn9Do1SXfrrTBoEGzevONYkyZ23CVYubLZC4r25eGJbRj0B9h777ADiwbvuskQw4dD8+Y+lziVBg60lcd77mnVEvfc0773olpJEJTN1no5XMm/2GXnUm65Bd8FLUG8RZ8BliyBCROsaNkuu4QdTXYZONATe0oEZbMn3l/ItLF9uPc2aLnMd0FLFG/RZ4A77oBGjaymjXNRtbV7L66ZdRpdu1p5D98FLXG8RZ/mCgttV53Bg8ErQ7go+8c/YOVKmDrVNrf/3y5oZS16n11Wa57o09xjj0FxMVxxRdiROJc8ixbB3/9u60SOOSY46LugJYwvmEpjpaU266BzZ5g2LexonEuOkhLo3RtWrICFC72mTW35gqkMNXWqrcAcNCjsSJxLnhEjbDHgv//tST5ZPNGnsZEjoWVLOPXUsCNxLjlWrIC//hV++UurUOmSwxN9mlq7FiZNgt/9Dho2DDsa5xKvpAR+/3vbAvPBB22tgksOH4xNUxMn2iDsueeGHYlzyTF8OMyYAaNHQ7uY+865RPEWfZp6/nnb7f7gg8OOxLnEmzULbrwRzjrLrlpdcnmiT0Pr19ussl//2i9nXfSsXw/nnGPF4R56yP+Pp4J33aShF1+0zb9POy3sSJxLLFW4+GL48kuYORN23TXsiLKDJ/o0NGGC9Vn26BF2JM4l1siRtqfC8OHQs2fY0WQP77pJMyUl8Oab0L+/7VPqXFTk58Pll0O/fnDttWFHk108laSZ+fOtD/Ooo8KOxCVdFpXgLSy0MacOHeC///VGTKp5102amTHD7j3RR1x+9pTg3brVkvyGDTBlCrRoEXZE2cc/V9PM9Om2KbVvVxdxWVKCVxX++EcrcfDEE9C1a9gRZadqW/Qi0gF4AtgdKAVGquoIEXka2D94WnNgvap2i3H+J8AGoAQorqzoTtbLz0ffzGPGtD8z4OT6YUfjki1LSvDee69VYL3hBmvVu3DE03VTDFyjqu+LSDNgjohMUdUzy54gIncD31XxM/qo6ro6xhpdwWX8qqKWrCsdymGtPgb2Cjsql0xZUIJ34kS46io4+WS46aawo8lu1SZ6VV0FrAq+3iAii4F2wCIAERHgN0DfJMYZbcFl/MLSLgAcuD4fT/RZIMIbvM+aZUXKevTwwdd0UKN/fhHpBHQHZpc7/HNgtaour+Q0BV4XkTkiUmnBXREZJCIFIlKwdu3amoSV+YLL+IXyMwB+eup+4cbjXCxxzhJautSqUbZvb4v/mjRJUXyuUnHPuhGRpsBzwFWq+n25h84Gnqri1N6qWigibYApIrJEVWdUfJKqjgRGgm08Em9ckRBcxi+6shGtlm+nzYDDwo7IuR+Kc5bQV1/ZGpCcHHj1Vd/+Ml3E1aIXkfpYkh+rqs+XO74T8Gvg6crOVdXC4H4NMAHw9Z6x9OrFogbdOfAgH4h1aSiOWULffQcDBsCaNTB5su2O5tJDtYk+6IMfBSxW1XsqPHwssERVv6jk3J2DAVxEZGfgeGBB3UKOrsJCn1bp0lTZLKGcnJizhDZuhBNPhHnz4Jln4DC/KE0r8XTd9AbOA+aLyIfBsb+o6svAWVTothGRPYBHVfVEoC0wwT4r2An4r6q+mqjgo2bNGr/UdWmqillCmzdbn/zs2fD005bwXXqJZ9bNTCBmIVFV/V2MY4XAicHXHwNeUT0OmzfDpk3Qpk3YkThXiRizhLZuta0up0+HJ5/0iqvpyksgpImyiUae6F2m2LYNzjgDXn/ddok655ywI3KV8dmtaWJdsJysZctw43AuHkVFcOaZ8NJLtt/rBReEHZGrirfo00Rpqd3v5O+IS3ObNll3zZQpVuJg8OCwI3LV8bSSZjS7VhC4DPPdd3DSSbbydfRob8lnCk/0acL3zXTpbt06Www1dy6MG2f98y4zeKJ3zlVr1So49lj4+GOYNMmnUGYaT/RpokEDu9+yJdw4nKto0SI44QT4+mt4+WXo0yfsiFxN+aybNNG2rd2vXh1uHM6Vl5cHRxxhUymnT/ckn6k80aeJVq2slKsnepcuxo6F44+HPfaAd96BQw8NOyJXW57o00ROji2W+uqrsCNx2U4Vbr0Vzj0XeveGt9+GPfcMOypXF95Hn0Z2390KmzkXlqIi2+N19GgYOBBGjYKGDcOOytWVt+hTqZqNG/bfHxYuTHFMzgUKC61e2ejRtsfrf/7jST4qvEWfKnFs3NC9u1X/++YbaNEipDhdVsrPt827N2ywMsOnnx52RC6RvEWfKrE2bqjQwu/e3Z46d25oUbp0EOeWfYny6KNw1FG25d8773iSjyJv0adK2cYNZS36li1/1MLv1s1a+B984NPYslacW/YlQlERXHUVPPSQza556im/kowqb9GnStnGDbfcYvdff/2jFn6bNtCpk81Xdlkqji37EmHFCpsf/9BD8H//ZwuhPMlHl7foU6nixg3lW/jB1mwDBsBjj9kK2caNsRZejF19XERVvPKrsGVfIowbB4MGWaXUCRPglFMS/itcmvFEH5ZKtmY76SS4/347fELz1F3GuzRRxZZ9dZKfz+bXZ3LV+7/lkRfa0quXJXzfozg7eKIPSyUt9aOPtkGxyZPhhHZ5P76M90QffTG27KuT/HwW9bmUM4vGsIC2XHfelwwb1Y769RP3K1x68z76MJQNuN1wg92Xm13RqBH06wfPPgvbevexlnxOTtIu4120lZbCiFs3cmjR26ymLa/WO5HbD3jCk3yW8UQfhmoG3C6+2GreTFh1+A8HcL01nxopnt6YLJ98Yu2IqyYfxzH18phb7xD6NczzBkMWqrbrRkQ6AE8AuwOlwEhVHSEiNwEXA8G21vxFVV+OcX5/YASQAzyqqsMTFHvmqmbArV8/6NzZ9uI8My/Bl/Guaimc3pgsqla64OqrbUObUaPggi7NkemX+qB+loqnj74YuEZV3xeRZsAcEZkSPPZPVb2rshNFJAe4HzgO+AJ4T0ReUNVFdQ08o1Uz4FavHvzhD3DddVYS4ac/DSVKk22zfmJdbWXQ6161yq4IJ0+2tRiPPVZWkKwXHJE5r8MlVrVdN6q6SlXfD77eACwG2sX583sAK1T1Y1XdBowDTq5tsJHSqxcMHVppErnwQuuvv/32FMdVXhVjCZGUnw+ffWbzDjNsXKS01ObEH3AATJsG//43TJ3qVSedqVEfvYh0AroDs4NDl4nIPBEZLSK7xTilHfB5ue+/oJIPCREZJCIFIlKwdu3aWE/JKq1b26rFsWNtpWwoUrR4Jy2Ufag98oj1fVx8ccZ028yfD0ceCZdcYjXj586Fyy+3K0PnoAaJXkSaAs8BV6nq98CDwN5AN2AVcHes02Ic01g/X1VHqmququa2bt063rAibcgQ2G0368IJRdlYQoa1bmul/IdaSYlNME/zJL95s10UHnIILF8OTzxhrfh99w07Mpdu4kr0IlIfS/JjVfV5AFVdraolqloKPIJ101T0BdCh3PftAa+4HqfmzeH66+H11+2WchXLNqR54quTDPtQe+UV6NoVhnnKxP8AAAtgSURBVA+H886DxYvtXmI1rVzWE9WYDewdTxARYAzwjapeVe74T1R1VfD11UBPVT2rwrk7AcuAY4AvgfeAc1S1yqrrubm5WlBQUIuXEz1bt9ofdEkJzJsHzZqFHVGEZcDA85Il8Kc/WaLfbz94+OG0/0xyKSIic1Q1N9Zj8bToewPnAX1F5MPgdiJwp4jMF5F5QB/g6uCX7SEiLwOoajFwGfAaNog7vrok736oUSMYMwY+/RSuuSbsaCKumgHyMH37rY3Z/OxntrXfXXdZ33zaJfmIrEGImmpb9GHwFv2PDRkCd95p0+ZOPDHsaFyqFBfb+PANN9iGNBdfbD1pbdqEHVkMEViDkMnq2qJ3aWDYMOvCufBCmwHook0VJk6Ebt1sD9euXeH9962rJi2TPGTXLK0M44k+QzRsaNsMbtliFS6//z7siFyyvPEGHH44nHoqbN9udY/efNOSflrLsAHtbOKJPpkS3F954IH2R79oEZx9tl3Wu+h45x3r+Tj2WFvhOmqUrYw+7bQMmU2TTbO0MoyXKU6WJPVXHncc3HefLY657DJ44AFfGJPpCgosN77wgi2UGzHCNgZp1CjsyGoh0SWWXUJ4ok+WJNZMGTzYZuEMH26X9iNH2tWyyxyqMGMG3HabrZFo3hz+/ne48kpo2jTs6FzUeKJPliRvCXfbbfZjhw2zfvsxY/Aa4xlA1WZO3XabXfS1bQt33GEf3rvsEnZ0Lqo80SdLsraEC4jAzTfbvrJDh8KmTfDkk76gKl0VFcH48Tb/fd48Kzb2wAPwu98FewM7l0Se6JMpEf2V1azWvO46u9S/8kro2dM2e95//7r9Spc4q1ZZVcmHHoI1a2xA/Ykn4Kyz/ArMpY4n+nQW54DuZZdZAjnzTDjsMEskp5wSQrzuf2bPtlLB48fbMM2AAXDFFTajJiNm0LhI8fka6awGC1D69oU5c6BLF5t/ffnlsGFDjCf6EvWk+e47W9B02GE2D/6ll+xDeNkyePFFmzHlSd6FwVv06ayGA7odO9pMjiFD4N57bbreww9D//7BE3yJesKVlsL06TB6tK1x2LrV6tHcdx/89rc+ZuLSg7fo01ktFqA0amTzsN9+G3beGU44wcrXrlqFL1FPoI8/trdln33saurFF+GCC+C992zjj0sv9STv0ocXNYuwoiK49VbrqalfH648/Uv+75nD2G37mh+26DOgPG86+PRTeOYZK0VR9t/zmGOs/tCpp/rsGReuqoqaeaLPNLVIyitWwI03wlNPwa5NixlyxFtcdm0Tmh7TM/XdORn2ofL559YlM368lSgAyM21ge8zzvA9WV36qCrRex99JqllUt5nH9t7dsgQuP76nRj6Uh9uz7euhktz5rJvklbwJir+VCopsYT+8su2sGnuXDvevbutRD7jDNhrr3BjdK6mPNFnkjqWVTjoIOtLLpv698ADMGL7YPrX24tL693P8fXzaJDMioNJLAtRF6tX22fO5Mnw6qtW9z0nxzbcvuMO65bxfVhdJvNEn0kSVFahZ09r4d91l9XJeejeo/nl18ezW/1iTn10J37zvQ0wJnxBT5LLQsSrsNBmypTdliyx461bWwnoAQPg+OOt/oxzUeB99JkmCX3c27ZZYa3x42HSJKt136KFJb1jjoE+faBDh+p/TlxS3Ee/ZYuVHJgzxwZQ33rLxizAassceSQcdZSFk5vrlUBd5vLBWBe3rVt3JP3XXoN16+x42TTCI4+0DTC6dEmvJfyq8NVXsHSp1esvS+wLF1pPEUCrVnDEEZbYjzrKXodX/XRR4Yne1UppKSxYANOm2Q5HeXk7drZq0MDKLhx8sC0Q6tzZZqB07GgJNRkrQDdtgi+/tK6XwkL46CNL7EuX2urT8rtutWplLfRDD91x69DBV6a66PJE7xKiuNiS6ty5P7x99dUPn9e4sSX8tm1h11133HbZxe7LXwmUJV5V2LzZkvWGDXb7/nu7rV5tiT3W9okdO1oRt/K3Ll2gffskJPUMmxrqskudEr2IdACeAHYHSoGRqjpCRP4B/BLYBnwEXKCq62Oc/wmwASgBiisLpDxP9Jnl669tMdFnn9nt00/ttnat1X/5/nu7/+67Hd0olWnY0D4QmjXbcd+mDbRrB3vsYbeyrzt0sNW/KZEBU0NddqvrPPpi4BpVfV9EmgFzRGQKMAUYqqrFInIHMBQYUsnP6KOq62oTvEt/LVva7ZBDqn5eWau9LNlXbGM0aZJe/f4/kKZTQ52LR7WJXlVXAauCrzeIyGKgnaq+Xu5p7wCnJydEFxUiKWyBJ1qaTA11rjZqNI9eRDoB3YHZFR66EHi6ktMUeF1EFHhYVUdW8rMHAYMAOnbsWJOwnEu+JO8Y5lwyxT0YKyJNgenArar6fLnj1wO5wK81xg8TkT1UtVBE2mDdPZer6oyqfpf30TvnXM1U1Ucf1/IQEakPPAeMrZDkzwdOAgbGSvIAqloY3K8BJgA9aha+c865uqg20YuIAKOAxap6T7nj/bHB11+p6uZKzt05GMBFRHYGjgcWJCJw55xz8YmnRd8bOA/oKyIfBrcTgfuAZsCU4NhDYF01IvJycG5bYKaIzAXeBSar6quJfxnOOecqE8+sm5lArKUnL8c4VtZVc2Lw9cfAwXUJ0DnnXN14CSfnnIs4T/TOORdxnuidcy7iPNE751zEeaJ3zrmI80TvnHMR54neOecizhO9c85FnCd6Fw35+XD77XbvnPuBGpUpdi4t+e5PzlXJW/Qu88Xa/ck59z+e6F3mK9v9KSfHd39yLgbvunGZz3d/cq5KnuhdNPTq5QneuUp4141zzkWcJ3rnnIs4T/TOORdxnuidcy7iPNE751zEeaJ3zrmIE1UNO4YfEZG1wKdhx1EDrYB1YQeRQtn2eiH7XrO/3syzp6q2jvVAWib6TCMiBaqaG3YcqZJtrxey7zX7640W77pxzrmI80TvnHMR54k+MUaGHUCKZdvrhex7zf56I8T76J1zLuK8Re+ccxHnid455yLOE30NichoEVkjIgvKHWshIlNEZHlwv1uYMSZSJa/3JhH5UkQ+DG4nhhljIolIBxF5U0QWi8hCEbkyOB7J97iK1xvl97iRiLwrInOD13xzcLyziMwO3uOnRaRB2LEmiif6mnsc6F/h2HXAG6q6L/BG8H1UPM6PXy/AP1W1W3B7OcUxJVMxcI2qHgAcDlwqIgcS3fe4stcL0X2Pi4C+qnow0A3oLyKHA3dgr3lf4Fvg9yHGmFCe6GtIVWcA31Q4fDIwJvh6DHBKSoNKokpeb2Sp6ipVfT/4egOwGGhHRN/jKl5vZKnZGHxbP7gp0Bd4NjgemfcYPNEnSltVXQX2hwO0CTmeVLhMROYFXTuR6MaoSEQ6Ad2B2WTBe1zh9UKE32MRyRGRD4E1wBTgI2C9qhYHT/mCCH3geaJ3tfEgsDd22bsKuDvccBJPRJoCzwFXqer3YceTbDFeb6TfY1UtUdVuQHugB3BArKelNqrk8USfGKtF5CcAwf2akONJKlVdHfyhlAKPYH8okSEi9bGkN1ZVnw8OR/Y9jvV6o/4el1HV9UAeNj7RXETK9tFuDxSGFVeieaJPjBeA84OvzwcmhRhL0pUlvMCpwILKnptpRESAUcBiVb2n3EORfI8re70Rf49bi0jz4OvGwLHY2MSbwOnB0yLzHoOvjK0xEXkKOBora7oauBGYCIwHOgKfAWeoaiQGMCt5vUdjl/QKfAL8oaz/OtOJyJHAW8B8oDQ4/Bes3zpy73EVr/dsovseH4QNtuZgjd3xqjpMRPYCxgEtgA+Ac1W1KLxIE8cTvXPORZx33TjnXMR5onfOuYjzRO+ccxHnid455yLOE71zzkWcJ3rnnIs4T/TOORdx/w88Vux9HlF44AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "t = np.linspace(0, 2 * np.pi, 100)\n",
    "p = ransac_model.params\n",
    "xe = p[0] + p[2]*np.cos(p[4])*np.cos(t) - p[3]*np.sin(p[4])*np.sin(t)\n",
    "ye = p[1] + p[2]*np.sin(p[4])*np.cos(t) + p[3]*np.cos(p[4])*np.sin(t)\n",
    "plt.clf()\n",
    "plt.plot(data[inliers,0],data[inliers,1],'r.')\n",
    "plt.plot(xe,ye,'b-')\n",
    "plt.plot(p[0],p[1],'bo')\n",
    "plt.axis('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we first determine position and radius of a reference sphere using laser scanner data.([SphRANSAC](https://nbviewer.jupyter.org/github/gyulat/Korszeru_matek/blob/master/SphRANSAC_en.ipynb)). Then we solve a more complicated geometrical problem: determine parameters of a best fitting right circular cylinder using a point cloud with outliers ([CylRANSAC](https://nbviewer.jupyter.org/github/gyulat/Korszeru_matek/blob/master/CylRANSAC_en.ipynb))."
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}