single_curve_reconstruction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "illegal-diana",
   "metadata": {},
   "source": [
    "# Single Curve Reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "indian-blind",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.insert(0, \"../\")\n",
    "\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from skimage.morphology import dilation, disk\n",
    "\n",
    "from curve_reconstruction.linear_model import LinearModel\n",
    "from curve_reconstruction.irls import IRLS"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "later-reset",
   "metadata": {},
   "source": [
    "### Generate line with outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "convertible-gasoline",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "w,h = 100,100          # image width and height\n",
    "sigma = 1              # gaussian noise scale\n",
    "a1_0, a1_1 = 30, 0.5   # line intercept and slope\n",
    "\n",
    "# Generate line\n",
    "A = [[a1_0, a1_1]]\n",
    "model_init = LinearModel('line', A)\n",
    "x = np.arange(w)\n",
    "y = model_init.predict(x,scale=sigma)\n",
    "y = list(map(int,y))\n",
    "\n",
    "# Generate image\n",
    "im = np.zeros((h,w), dtype=bool)\n",
    "im[y,x] = 1\n",
    "# add outlier\n",
    "im[10:15,80:85] = 1\n",
    "im = dilation(im,disk(1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dedicated-express",
   "metadata": {},
   "source": [
    "### Ordinary Least Square Fitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "overhead-floor",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Model data\n",
    "y_im, x_im = im.nonzero()\n",
    "X_im = model_init.basis_functions(x_im)\n",
    "# OLS Parameters estimation\n",
    "a_ols = np.dot(np.dot(np.linalg.inv(np.dot(X_im.T,X_im)), X_im.T), y_im)\n",
    "# Curve reconstruction\n",
    "model_ols = LinearModel('line', [a_ols])\n",
    "y_ols = model_ols.predict(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "knowing-america",
   "metadata": {},
   "source": [
    "#### Plot Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "caring-amendment",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f91f07334d0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcEAAAHWCAYAAAAPaDLLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAnAklEQVR4nO3deXhcV53m8fcny/Ie77vjTTKBkKSB2CFAFidkkSGkYRiGZeBJyEBDDA09DNBPeoCYMNPrdLNMj92d0HSagWkGCDMwgG1iYpOEEJIYTFiDJe92HMvyvtvSmT/OLdWiulelUl3Vcr6f59GjpOrce08trlfnnvO7Zc45AQAQoqZqdwAAgGohBAEAwSIEAQDBIgQBAMEiBAEAwSIEAQDBIgRRMjNbaGbOzB4suP3B6PaF1ekZAJSHEKxjZrbUzP7ZzLaZ2WkzO2ZmvzSzvzGzudXuXyMxs+VR0G+qdl8Gw8xWRf1eXsa2o83so2b2UzM7ambnzOx5M9tsZn9vZtdXvsfA8GqudgcweGZmkv5S0sclXZD0sKRvSGqR9GpJH5W00szucM59cxi6dE/Un73DcCwMAzMbL+lHkl4hab+kh6Lf4yX9gaQ/kjQpagPULUKwPn1SPgB3SLrNOffr3DvN7M2SviLpa2Z2s3NuY5qdcc49L+n5NI+BYfcn8gH4A0lvcM6dy73TzCZLekkV+gVUFKdD60w07/ZJSecl3V4YgJLknHtI0n+UNELSGjNrytn+zuj02J1m1m5mm6JTXS6nzQQz+zsz22NmZ8zsd2b2EcW8X4rNCebOH0b//TUzOxjt7xkzu63Ifiaa2cfM7JHo2OfMrMvMvmNmr4o5tosewywz+6KZ7TWznujx/Wt0f9HTdmb25uj+vy/6ZA+BmY01s3vMbIuZnTSzE2b2EzN7e5G2LWb2QTP7vpntNLOzZnbIzDaY2YqY/V8RPb4dUfsuM/uZmX3OzEZGbXZIujfaZGP0WF3ua53g1dHvNYUBKEnOucPOuSeK9Gummf2Tmb0QnaLfYmZ35JxOXlXQfkfUz2KPseipXDN7o5l9xcx+Hz23J6NTtB/Kfa/ntM+8Pxeb2R+b2bNR3zbltJliZn9hZr+N7jtqZj80s1uK7K8lOtbPzOywmZ2KHse3zeymYo8FtYuRYP15t/zr9nXn3C8T2n1R0qckXSLpekmFo8F/K6ld0lpJ/yBpgSSZ2ShJP5S0TNIvJH1V/rTXJ6P9DNYCSU9J2ibpf0qaIumtkr5tZjcVjFJfIum/SnpU0vckHZY0X9LtklaY2Rucc+uKHGOKpCclnZD0LUm9kl6QtEbS2+RP3RU7bfe+6Pc/lPG4YpnZJEmPSHq5pJ9J+pL8HxC3SvpfZvZS59wnCvr/eUlPyJ/a7pI0W9IbJH3fzN7rnPtizv6vkPRTSU7SdyRtl3SRpDZJKyV9Qv6PpM9JeqP86/Yv8mcOStUd/X5RqRuY2bToMSyW9Hj0M1v++f3BII49kL+Uf41/Kn8KfqKkG+Wfw2WS3hWz3eclXSv/3vq+pJ6o3wskbZK0UNJjktZJGifpNknrzOx9zrkHcvbzoKS3S/qVpC9LOi1pjqRr5P9NbajEg8Qwcc7xU0c/8gHlJL23hLZfjdp+Iue2O6PbeiW1F9nmz6L7H5LUlHP7IkmHovseLNjmwej2hTm3LYxuc5LuLWh/a3T79wtunyhpWpE+zZO0T9Jvi9yXOcaXJTUXuf9Xks5Imlpw++LoOfhxic/78ug4m0pom3k+Pl5w+2j5D9heSS/LuX2UpHlF9jMx6v8hSWNybv/baP9/WGSbyQWv26qo7fJBvs9ui7Y7K2m1pNdLmj3ANvdH23y24Pal8qHsJK0quG+HpB0x+yvad0mtRdo2yQe9k/TKmNdjr6RFRbbdFL0mbyu4fZKkLfIhNzPnNemV9IykEUX2NbXYY+Gndn84HVp/Zke/d5fQNtNmTpH7vu2Kj6reLf+P/OPOud7Mjc657ZK+MJiORnZK+i+5Nzjn1kvaJemqgtuPOucOFu7AObdH0jclvdjM5hc5xjlJH3XOXShy3xr5kLmz4Pb3SjJJ/1jawyiNmU2V9E5Jzzjn/jr3PufcGUl/Gh33HTm3n40eowraH5UfRU6WH+EUOl1km8O5r1u5nHPflfTh6Bh3S/qupH3mV4d+1cyuy20fnYL995KOy4dX7r6ekf+DrCKcc51FbuuVH+lJ/o+sYv46eh/3MbM/kB8pP+Sc+1rBPo/In04eLenNmZvlX7+z8v9OCvvRXXgbahunQ8P1VOENZjZB/pTa7mIfNPJ/Md9b5PYkW5xzPUVu3y2p3zyfmb1G/sP3VZJmyK94zTVXPkBz7XDOHYg5/pflT5/9kfwIKvOBfaf86davl/QoSrdMfi623/xXZGT0O29RiZm9VNLHJF0n/4fO6ILtckte/rf8c/R/zeyb8qfffhzzmpXNOfcFM/uipJvl5whfHv1+h6R3mNlnnHOfipq/WNJYSY9F4V1ok6Q7KtGv6A+Nj0l6nfyIflxBk7jyoH7veWXfgxNjXq/p0e+XSJJz7piZ/T/5U9VbzOwh+VOoP3XOnSr5QaBmEIL1Z7/8P8iLS2ibabMvZj+FJka/X0g49mAdibn9ggoW2pjZm+RHfGfk58Y6JZ2U/4t7ufxf7KMG0y/n3HEz+4qk95vZDc7PQd4uaZakz0Wjs0qaGv1epuKjt4zxmf8ws6vl5xCb5U93f0fSMUWnTSX9oXIet3PuKTO7VtJ/lp/bfVe0n+ckfdo5968VeiyKPti/Hf3IzFrkR9Gfl/RJM/uWc26L0nnv9BPNtz4tf3r+Kfk/cg7Jv58myf9xUOw9EteHzOt1c/QTZ3zOf79VfkT/Dkmfjm47E/1B8lHnXNxzgBpECNafxyXdIOkmSQ/ENTKzEfLBIUk/LtKk2ArBzF/wM2N2O6u0LpbtM/KnNpc6536be4eZ/aPiF+YMtNpxjaT3yy+E2ajsgpj7y+9qrMxz+Fnn3EdK3OYTksZIusE5tyn3DjO7Rz4E8zjnfiLptmgh05XyCzL+WH7hTZdzLpXFGc6vFP0fUXC/U35ByhaV/97pVf/RfsakIre9Rz4AP+2cW5V7h/kVxB+O2ZeU/J7/sHOupNP9zrnT8qd8V5nZxfKj9zvln4+F8otvUCeYE6w/D8qvantTdAotzl3yc4HPqcSCZufccUkdkuaaWWuRJssH1dPBa5P0myIB2CS/8q4szrln5f8QeJOZvVL+D4hHC49TIU/Jf7AP5oOwTdKhwgCMJK7IjeYTn4hOS34oujk3NDOnokcMoj+lOB79tuj37ySdkvQyM5tYpP3ymP0cljQzU9ZRYGmR29qi3w8Vua+c1ctPRr/LCi7n3G7n3Ffl5yE7JF0Tna5FnSAE64xzbpukP5efW/qOmV1a2MbM3ih/uqpH0t2DXCjxz/Lvi7+y/PrCRcp+yKZlh6QlZta3kMfMTP6v7n6Pc5DWyI84HpL/4K5oWURGNDf5VUlLzeyT0Yg8j5m1Rs9nxg5JU6LSh9x2/0FFFnmY2avNbEyRw2dGYblzU5mFGsUWFMUys/dHo71i971Y0lui/31Ukpxz5+Uf9wQVLIwxs6Xyi2aKeUr+jNS7C7a5U9JrirTfEf1eXtD+5fJXLhqUaNHOY5L+jZndVayNmV1uZjOi/55uZpcXaTZO/pTpBfmzGagTnA6tT6vk/9F9RNIvzGy9pF/LB+OrJb1SflXf293grxbzt/K1ZW+W9LNo35Mk/Tv5D7zbh979WJ+VD6efRwsOzst/EF4qKbMYoVzfiPY/V9JB+XrCcrzYCi4gnmNXNCL7oKQlku6T9C4ze1x+rmyO/HzuMvk6s8xKxc/Jh93jZvZ1+VN0S+VHv9+Un/fL9XFJN5rZY9E+Tkh6qaQV8iOr3NO8G+VHpn9hZpdF98s5l7dit4h2+Qst7JAfRe+Wn2tbEvV1pKQvOOeeztnmzyS9VtKfRMGXqRN8q3xdXrH3zn+XD8A1Zvba6Dgvk1+w8l35Uo1cX5ZfFPM5M7tB0taoT7fJv6ZvHeBxFfMO+TnZfzKzD8nXHx6RL825QtJlUX8OyL9/fm5mv5T0bNTfi6Ljz4qek+OFB0ANq3aNBj/l/8iXGPyL/AfhafkPw19J+m8qUncWbXOn/NzInQn7vUjS38nXVZ2RP9X1n+RX4g22TvDBmGNs8m+/ov3bIr8g5qCk/yPpcsXXjJVUuxe1/WzU/m/KeK6XK1uTGPezJad9i3wYPiEfamflV7X+UP6SZIV1i7fJn5o7Lv8B/ANl55ryXi9Jt8iP2H8T7fuk/GnvL0haUKTv71S23s0Ve96LbPOi6DVfK3+a72TOY/iW/OX6im03S76soys63pboMWSev1VFtrlG/g+sU/ILgr4nHz5xr/ml8ouHDkT92iw/V1j0Paci788ifZggH+Kb5f8dnZb/d/U9+ZXF46J2k+QvQvGI/L+Ps/KXDNwk/4eNVftzgZ/B/Vj0wgINL7pM1nWSLnHOba1yd4Ji/tJnG1VkQQtQTcwJIghmdpX8won1BCCADOYE0dDM7G75eZzMlXAGW+wPoIERgmh0fyq/wGGbpHc554pdNQRAoJgTBAAEizlBAECwEk+HWmlfvgkAQM1yzlncfYwEAQDBIgQBAMEiBAEAwSIEAQDBIgQBAMEiBAEAwRryFWOuvvpqrVy5UlOnTlVTU5iZ2tvbq+7ubq1evVpPPvnkwBsAAGpC4hVjBqoTvPrqq/WpT31K06dPr3jH6lFXV5fuu+8+ghAAakhqdYIrV64kAHNMnz5dK1eurHY3AAAlGlIITp06tVL9aBg8JwBQP4YUgqHOASbhOQGA+sEnNgAgWIQgACBYhKCkBx54QK973et01VVX6fbbb9eyZcvU0dEhSTp//rzuv/9+Pffcc1XuJQCg0oL/Zvnf/OY3uv/++7Vy5UpdeeWVGj16tM6ePat58+ZJ8iH4wAMPaPbs2brkkkuq3FsAQCUFH4I7duyQJL3lLW/R+PHjq9sZAMCwCvp06KpVq3TvvfdKkm644QYtW7ZMmzdvzjsdev3110uS7rvvPi1btkzLli3Tvn37qtZnAEDlBD0SfM973qOZM2fqS1/6ktasWaNRo0bp0KFDeW3WrFmju+++W3fddZeuueYaSdK0adOq0V0AQIVVPASXfW9ZpXc5aE+//umS2s2bN69v7u/SSy/V2LFjtXnz5rw2l156aV/byy+/vLIdBQBUVdCnQwEAYSMEAQDBqvjp0FJPRQIAUG2MBAcwcuRISdK5c+eq3BMAQKUFvTq0FCNHjtScOXO0YcMGtba2qqWlRUuWLOkLRwBA/WIkWIJ77rlHR44c0Qc+8AHdcccd6urqqnaXAAAVMKRvll+7di01cwUOHjyoFStWVLsbAIBIat8sDwBAPSMEAQDBIgQBAMEiBAEAwRpSCPb29laqHw2D5wQA6seQQrC7u7tS/WgYPCcAUD+GFIKrV6+mZi5HV1eXVq9eXe1uAABKNKQ6QUm6+uqrtXLlSk2dOlVNTWFOMfb29qq7u1urV6/Wk08+We3uAAByJNUJDjkEAQCoZRTLAwBQBCEIAAgWIQgACBZfpQQANSBpfcZQmMVOh0GMBAEAASMEAQDBIgQBAMFiThAAqiCtOcCBjsMcYT5GggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGA1V7sDABAiM8v7f+fcsBwH+RgJAgCCRQgCAIJFCAIAgkUIAgCCxcIYAKgBLGCpDkaCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgNQ+msXOu321mVrHOAAAwnBgJAgCCRQgCAIJFCAIAgkUIAgCClbgwpthCmIHasFAGAFAvGAkCAIJFCAIAgkUIAgCCRQgCAIJFCAIAgkUIAgCCRQgCAII1qAtol4K6QQBAvWAkCAAIFiEIAAgWIQgACBYhCAAIVuLCmFIWtQx0ke1SLsJdznFYcAMAGCpGggCAYBGCAIBgEYIAgGANuVi+cG6ulDnAQmltw7whACAJI0EAQLAIQQBAsAhBAECwCEEAQLAq/i0SlVgoU2xBSyUWz7BQBgCQi5EgACBYhCAAIFiEIAAgWBWfEyxUzvxesfuHqyifeUMACAcjQQBAsAhBAECwBg7BeZLGpN8RAACGW/KcoEl6m6RxkvZK6oh+9koa/JRcdrdlzO+VU/NHbSEAIIklBYXNMaf3FbnjlKRtyobiiaF1olIF9ZXYbznHAQDULudc7Ad5cgguMKebJc1V8onT/ZK2ygfibkm9g+7g4DYQIQgAKE35IWjm7xwjabGktuhnQsLRzip/lHi0pA4O3Kh/31LZbznHAQDUrqGHYKGZygbifEkjEo7epWwg7pR0oWgHE3YQ27cB2xCCAIDKh2CuUZIWSWqVtETSpIS25yXtUDYUu4s3K6VYvlLSCmAAQG1INwQLTZUfIS6RtEDSyIS2h5UNxO2SzvmbCUEAQKUMbwjmapYPwkwoTkto2yNpl6St0rPfelaXzbgsL2wIQQBAOaoXgoUmKTuXuEj+VGqMORPm6NbWW7WibYVuWnyTpoydUtGuZBCCANDYaicEc42QdLGyoTgroW2vpD3Knjp9XmUV66d1YW5CEQBqV22GYKEJ8otr2qLfSZdqOympUz4QO6P/LwEhCADhqY8QzDuwfIF+ZpQ4N7otzj5lR4l7FFusTwgCQHjqLwQLjVV+sf74hLZnlF+sfyx7FyEIAOGp/xDMZcoW6y+Rn1dMuKTbS6e/VO1t7Wpva9e186/VqOaE1TiqTCiWguAEgOHRWCFYKFOsv0Q+GCfGNx07cqxuXHSj2lt9KLZOae3XhhAEgMbS2CFYaLr6Tpu2vKhF53rOxTZtm9LWF4jLFy7XuJZxhCAANJiwQjDHibMntGnHJq3vXK91Heu09dDW2LYtI1p03YLrtGHNBj+X2OVvJwQBoL4FG4KFj63zUKfWdazT2o612rhjo06dPxW/8VFJHdJDf/WQXrvotZo4OuE8axkIQQAYHsGGYKHcx3r2wlk9vutx3fT+m/zp05kJG/bKf09itOK0Z2+PmizpCxbTUcqpWsIVAPIRgpHE0LhI+cX6o+P3M2PcDN3aeqva29p18+KbNX3c9BR62x8hCACDRwhGSg6NJvkC/SXygTg3fp8m09I5S/vKMK6ae5Wam5or1OPkvhKCADAwQjBSdmiMky/Wz4TiuPimk0dP1s2tN6u9tV23tt2qORPmlNnb/ghBABg8QrCSTNJsZa9eM0+JxfpXzLyirwzjNfNfo5YRLcPSzQxCEUDoCME0jZa+sfkbWrt1rdZ3rtfe43tjm45vGZ9XrL9o8qLUu0cIAggdIZiyzHPonNOvu36ttVvXal3nOj228zGd7z0fu92Lpr5I7a3tWrFkha5fcL3GjEz66ozyEIIAQkcIpizuOTxx7oQ2bt/YV5u4/cj22H2Mbh6t6xdc37fA5pKpl1QkwAhBAKEjBGvFFPXNJY65dIxOXzgd23TBxAV9gXjjoht10aiLyjokIQggdIRgDTp17pQe2/VY3yjxdwd/F9u2ualZ18y/pm8u8YqZV5QcboQggNARgjWo8HnfeWRn3zVON2zboOPnjsduO2v8LD9KbG3Xza03a8qYKbFtCUEAoSMEa1DS836+57ye2P1EXyj+fP/PY9s2WZOumntV3yhx6ZylGtE0ou9+QhBA6AjBGjSYb6fYf2K/1nes1/pO/3Po9KHYtlPHTNUtrbeova1dt7TeotkTZleiuwBQtwjBOlXstenp7dHm5zf3lWE8tfcp9bre2H28fNbL+xbYvGreqzRyxMgBj8voEUAjIQTrVCmjxe5T3Xp428N9p073n9gf23ZCywTdtPimvlCcP3F+0XaEIIBGQgjWqcF+oa9zTs++8KzWdvir1zy+63Fd6L0Q2/4l017SF4jXLbhOo5v9V2cQggAaCSFYp4b6rfbHzh7TI9sf0bqOdVrXsU47j+6MbTumeYyWL1yu9rZ2ffj1H5a6h3RoAKgZhGCDKHytio3Y4l5P55ye636ury7xRzt+pLM9Z2OPtWjSIrW3tWtF2wrdsOgGjW8Zn3hsvtECQK0iBBvEUEKw0Knzp/Tozkf7Ftj8vvv3sW1HNo3UtQuu7SvDuGzGZWpqyv/qDEIQQK0iBBtEJUOw0LbD27S+Y73WdqzVI9sf0cnzJ2PbzpkwR/t+tE/qkLRN0hlCEEDtIgQbRJohmOtczzn9eNeP+xbYPPvCs/GNeyXtkT5zx2fU3tauV8x+hZqsKbZ/ADDcCMGAVCIEC+09tlc/6PyB1nas1cPbHtaRM0di204bO023tt7aV6w/Y9yMvvsIRQDVQAgGJI0QzHWh94Ke3vt03wKbZ/Y9I6f4Y145+8q+MoxrF13rR44AMIwIwYCkHYKFuk526eFtD2tdxzqt71yvAycPxDc+I6lTfi6xU9Kx4ekjgLARggEZ7hDM1et6tWX/lr66xCd2P6Ee1xO/wQvygdghaZekhKYAUC5CMCDVDMFCR88c1YZtG3wodq7TnmN74hufk7Rd2VA8PDx9BND4CMGAlbKidKBtytXvWNMltUU/CyQ1J2zcrWwg7pB0viJdAhAgQjBgNRWCuUZK3/3Vd/tGiR2HOuLbXpC0U9lQ7KpI9wAEghAMWM2GYMFxOg51aH3Hen3w8x+UFkpqSdjwqLKBuE1S/NXfAIAQDFm9hGDeNiMkzZc/bbpE0ox+zbJ6Je2WD8St8otteNcCyEEIouKG6+o1ZiZdpOxc4mJJoxM2OKH8UeKpIXcBQJ0jBFFxwxqCuZokzVM2FOckbOwk7VU2FPeKUSIQIEIQFVe1ECw0TlKrfCC2Rv8f57Tyi/WPD7l7AOoAIYiKq5kQzGssabayo8R58iPHOPuVHSXuFsX6QIMiBJG64SrSH1Qojpa0WLrrz+/Sus512nd8X3zbs8ov1j8ylF4CqCWEIFJXkyEYcc7JOadfHvhl33cmbuzc6Fehxjmo/GL9C2V1F0ANIASRuloPwX77GWXSImVPnU5O2MF5ZYv1t8pfzQZA3SAEkbq6C8HC/UxVNhAXShqZsMPDyo4St8tf9xRAzSIEkbq6D8FczfLXNs2E4vSEnffIfwNGJhRfGHT3AKSMEERVpBGMVVmFOlH5xfqjEnZ0XPnF+qeH3DUAQ0QIoioaJgRzjVB+sf7shLa9yi/W3yeK9YEqIARRFQ0ZgoXGK79Yf2xC21PyRfpbo98nyzskgMEhBFEVQYRg3k4kzVV2lDg3ui3OPmVHiXvkR44AKo4QRE0YrqvMlKIioTeQMfJziEvkR4kTEtqekV9pmhklHk29d0AwCEHUhOBCMO+AkmYqO0q8WMnF+geUHSXuEsX6wBAQgqgJQYdgoVHKL9aflND2vPIv6XYo7c4BjYUQRE0gBBNMU36xfnNC20PKv6QbxfpAIkIQdSutRS+lBHLVjFR+sf60hLYXlF+sfyD13gF1hxBE3QoyBAtNVrYEY7GkloS2x5RfrH8m9d4BNY8QRN0iBAuMkDRf2VCcldC2V770IhOKz4tifQSJEETdIgQHMEE+DJfIjxLHJLQ9qfxi/VOp9w6oCYQg6hYhOAhNyi/Wn6P4Yn2n/GL9vaJYHw2LEETDGPQ3QsRsUwk1H5xjlT1t2iZpXELbM/Kjw47o97HUewcMG0IQDYMQLJPJzx/mFus3JbR/QfnF+j1pdxBIDyGIhkEIVsho+WL9JfKheFFC23PKL9Y/nHrvgIoiBNEwCMGUzFD2tOkCJRfrd8uH4VZJO+WvaAPUMEIQKDBc33BRl1rkr1qTOXU6JaHtBfkgzIwSu9LuHDB4hCBQgBAchCnKBuIi+SvaxDmibCBul3Q27c4BAyMEgQKEYJmalS3Wb5M/jRqnR9JuZUPxBVGsj6ogBIEChGCFTFS2WH+R/IKbOCeUDcROSadT7x0giRAE+iEEU9AkaZ78CHGJpNkJbZ18gX5usT6fNkgJIQgMINivcErTOGVHia3yxftxTitbrN8hP2oEKoQQBAZACKbM5C/jlplLnKvkYv398iUYHfLzilzSDUNACAIDIASH2Rj5C35nQnFCQtuz8l8LlZlLPJJ259BoCEFgAIRglc1UNhDny39lVJwuZU+b7pSvVQQSEIJABQxXUJajocK1RX6laSYUJye0PS9ph7Kh2J1251CPCEGgAgjBKpmqbCAuVHKx/mHlF+ufS7tzqAeEIFABhGANaJa/tmkmFKcntO2R/waM3GJ9BIkQBCqAEKxBk5S98PdiSaMS2h5TtgxjmyjWDwghCKSAUKwxI+SL9TNfDzUroW2v8ov194li/QZGCAIpIARr3ARlR4mt8mUZcU4qO0rsjP4fDYMQBFJACNYRky/Qzy3WT3qK9ik7StwjivXrHCEIpIAQrGNjlB0ltkkan9D2jLLF+h3yc4uoK4QgkAJCsEGY8ov1L1Zysf4B5Rfr96TdQQwVIQjUiLSCszD0ih2nnGAs3E8Q4TpKfqVp5uLfExPanlN+sf6htDuHchCCQI0gBOvQNGVXnC6Qr1WM061sIO6Qv6INqo4QBGoEIVjnRipbrL9E/mo2cS7Iny7NhGJX6r1DDEIQqBGEYIOZrOxc4iL5657GOar8S7qdSb13iBCCQIMrJ1xLCc7B7iNoI+S/ASMTijMT2vbKf09iJhT3i2L9FBGCQIMjBGvQRcov1h+d0PaksoHYKelU6r0LCiEINDhCsMY1yRfoL5EPxLkJbZ3yi/X3imL9ISIEgQZHCNaZcfJlGJlQHJfQ9rTyi/WPp967hkMIAg2mWoX6hF4KTP5i35kyjHnyI8c4LygbiLtEsX4JCEGgwRCCDWy0/Cgxs8DmooS255Q/SjySdufqEyEINBhCMCAzlA3EBUq+pNtB5V/SjWJ9SYQg0HAIwUC1SFqobLH+5IS2F5R/SbeDKfethhGCQIMhBCFJmqL8Yv2RCW2PKL9Y/2zanasdhCDQ4Eq5sksawVnKcQjOYdKsbLH+EknTE9r2qH+xfgMjBIEGRwiin4nKFusvVnKx/nH5Iv2t8gttTqfeu2FFCAINjhBEoib50otMGcbshLa98gX6ncoW69d5EhCCQIMjBDEo45V/SbexCW1PyY8Ot8oH44nUe1dxhCCAqi2mKcVAV68hSFNkkuYou8BmrpKL9Z9Xdi5xt+rikm6EIABCEKUZIz86zIwUJyS0Pav8Yv2jqfeuLIQgAEIQ5Zmp7ChxvpKL9buUX6x/IfXelYQQBEAIYuha5OsRM2UYkxLanld+sX53yn1LQAgC6KecoKmlICU4a8A05V/SLalY/7B8GG6VD8dzaXcuixAE0A8hiIpqVvaSbm3yARmnR/50aWaUeCDdrhGCAPohBJGqScq/pNuohLbHlA3EbZLOVLYrhCCAfghBDJsRki5WNhRnJbTtlbRH2VB8XkMu1icEAfRDCKJqJii/WH9MQtuTyl69pjP6/0EiBAGkotjnx0DhVE0EYw0y+QL93GL9pJdpn7KjxD0qqVifEASQCkIQFTdW+cX64xPanlF+sf6x4s0IQQCpIASRKpMv1l8iH4rzlXxJtwPKlmHskl+FKkIQQEoIQQyrUfJfC5U5dToxoe05+S8P/oXkfk0IAhgmtRJ65X6TBkFZR6Yrv1i/uUibRyX3w/gQLLYJAAC1ryv6+Yn81WoWKhuKU6M2Hcm7YCQIoKIYCaImTJGfR9wsuR5OhwIYJoQgak3SwhhOhwKoqDSK2MsJr3LDmKL7sCQtNgUAoKERggCAYBGCAIBgEYIAgGCxMAZAqiqxsKSai1NKWSjDYpr6xUgQABAsQhAAECxCEAAQLOYEATSE4bpSTSnHYR6xfjASBAAEixAEAASLEAQABIs5QQANoZQLaJcy71aJucVyLubNHGF1MBIEAASLEAQABIsQBAAEixAEAASLhTEAGlK5C0sG2q6UBTflLK5hIUx1MBIEAASLEAQABIsQBAAEizlBABiEUi6GXY7hmkekKD8fI0EAQLAIQQBAsAhBAECwmBMEgCGqRJ1gOdKaiwxpnpCRIAAgWIQgACBYhCAAIFiEIAAgWCyMAYAKq0QRezWFVFDPSBAAECxCEAAQLEIQABAs5gQBoAbU0jxiJb5YuF4wEgQABIsQBAAEixAEAASLEAQABIuFMQBQp0pZjFJO4Xs5C27qtcCekSAAIFiEIAAgWIQgACBYzAkCQAMrp/B9uNTCPCIjQQBAsAhBAECwCEEAQLCYEwSAgBWbh6vEPGElag2l9OcJGQkCAIJFCAIAgkUIAgCCRQgCAILFwhgAQJ7CxSilFLVXYjENxfIAAAwjQhAAECxCEAAQLOYEAQCJSpmrG2gesVLbDLSPwe6HkSAAIFiEIAAgWIQgACBYzAkCAIasUhfMHu59MBIEAASLEAQABIsQBAAEixAEAASLhTEAgCErpfB9oKL7SiyUGSxGggCAYBGCAIBgEYIAgGAxJwgAqLhyviC3nG1KmXtMmmtkJAgACBYhCAAIFiEIAAgWc4IAgLpVzjxiLkaCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgEYIAgGARggCAYBGCAIBgmXOu2n0AAKAqGAkCAIJFCAIAgkUIAgCCRQgCAIJFCAIAgkUIAgCC9f8BvUc+53yzfzEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax = plt.subplots(1,1, figsize=(8,8))\n",
    "ax.imshow(im,cmap='gray')\n",
    "ax.plot(x,y_ols,'-g',lw=3,label='fit')\n",
    "\n",
    "ax.set_aspect('equal')\n",
    "ax.axis('off')\n",
    "ax.set_title('Ordinary Least Squares', fontsize=20)\n",
    "ax.legend(loc='upper left', fontsize=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fixed-operations",
   "metadata": {},
   "source": [
    "### Iterative Reweighted Least Squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "diagnostic-ceramic",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initial parameters: horizontal line\n",
    "a0_0, a0_1 = 50, 0\n",
    "A_0 = [a0_0, a0_1]\n",
    "model_init = LinearModel('line', [A_0])\n",
    "y_init = model_init.predict(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "coordinate-screening",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select a noise model from Table 1 of Tarel et al. 2002, and scale suitable for the problem\n",
    "noise_model = 'cauchy'\n",
    "scale = 5\n",
    "A_irls = IRLS(X_im,y_im,A_0,noise_model,scale)\n",
    "# Curve reconstruction\n",
    "model_irls = LinearModel('line', [A_irls])\n",
    "y_irls = model_irls.predict(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "approved-thailand",
   "metadata": {},
   "source": [
    "#### Plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "decreased-keeping",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f9221a19cd0>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcEAAAHWCAYAAAAPaDLLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA9DklEQVR4nO3dd3wUdf7H8fc3ld67oAgqHaQEEEhiOe8sd6I/9dQr6ik2UMFCDwkJRUUUFUURz3qWu1PPO+489CwkoQsCFgRFjaCA9A4hJPP7Y3ZhZ7Ob7CabbJJ5PR+PPJLMZ+b7/e5see+0XWNZlgAAcKOYaA8AAIBoIQQBAK5FCAIAXIsQBAC4FiEIAHAtQhAA4FqEYBVljDnXGGMZYyZHeyw1nWc9L4xAOwuNMVX2miNjzGTPbT032mMBqopqH4KeJ7XlN629Z/qLURpWqarDGP1517XPT6ExZrfnxf9GY4yJ9hirm0gFcEXwCc3J0R5LOMrzZsQY08YYM8sYs84Yc9gYc8QYs8kYk22MmWaM6Rjp8SK64qI9AAS1QlIXSTujPZAAMj2/4yWdIekKSamS+km6M1qDKocukg5HexCILmNMd0nZkppI+lzSS5J2S2ohqb+kCZK+l/RttMaIyCMEqyjLsg5LWh/tcQRiWdZk3/+NMYMl5Ugabox5xLKs76MysDKyLKtKrmdUusdkB+Bky7Iy/YvGmA6SEip7UKhY1X53qD/Prhvvi/ANfrvvbvSb91fGmHeNMTuNMfnGmG+NMQ8bYxoFaDfP89PAGPOo5+8C764iz26UdGPMYmPMNmPMMWPMFmPMa8aYruGOMdAxQWPMek+7zYLc9rGeZe70m97WGPOkMeY7z+3cZYz5lzEmKcTVWiLLshbLDmwjqW+AcQ0wxrzps142G2PmGmPa+M33umf8Z/pNf8kz/UO/6fU990FOgD6vM8Z8bIzZa4w5aoz5yhiTZoxJDDBvwF2SxpjWxpgXjDHbPbvF1hhjbgh03/gtF2eMmWCM+cazvjcbYx4yxiT4zHOjzy67VL/HwGS/9kJafz7z9zXGLDDGHDDG7DfGfGCMOSfQvJEUzuMsnOeLzzKXGWM+NMZs9bS/xdi7KYd76u096zTV87/vOl0Ywk0Y5Pn9eKCiZVnfBXrDFGx9mwDHYE0ph0FMgF25xpgEY8ydxn6t+sFz23d7+rk4SDslvl555ulsjHnR83g6Zoz52bP+OwVor6UxZqYxZoMx5pDnebXBs3yHQGOoLmriluBCSY0kjZS0VtI7PrU13j+MMRmSJsve3fFvSdsl9ZR0v6RLjDHnWJa136/tBEkfyX63+L6k/ToZZimSxkn6WNJbkg5KOlPSVZIuM8YMtixrbThjDOAlSdMlXSdpdoD6DZKOSXrN53b28Yy1iaT3JL0tqZmkyyUtMsZcYVnWuyX0Ga4C33+MMTdJelZSvqR/Sdose70Mk/QbY8xAy7I2eWb/UNK1ki6Q9I1PMxd4fg8yxtSyLOuo5/9U2Y9h/3B8XtKfJP0o+77YK2mgpCmSLjDGXGhZ1vGSboQxpoWkpZJOk72Vu0RSK0lzZK/PkrwmKVnSf2U/Ri6RNEb2brU/eeZZI3u3coakHyS96LP8Qp9xhLP+ZIwZJOkD2Y/VtyVtlHS2p82PShl3mZXhcRbO80XGmFslzZW0TdJ82YcJWsh+zv5J9v2yV/Y6vVH2/ea7NZcXws3YJamtpLNkH44I5XZXxvpuIjuYl0j6n6QdklpL+o2kd40xt1iW9VyA5YK+XhljLvKMN172+two+7b/n6RLjTHnWZb1qWfeOpIWS+ro6X++7De8p0kaKulNSd9F6LZWPsuyqvWPJMu+GY5p7T3TXwyyzHme+hJJjfxqN3pqs/ym53mmfyCpboA2W0iqH2B6L9lP8P+GOcZzPfXJPtPaSiqUtDLA/Eme+d/ymRYn+8F9VFKq3/xtJP0kaaukxLKua8/0FM+48iW19pl+luxQ3ijpFL9lLvAs8w+faR08ffzdZ1onz7T3Pb8v8KnN8kxLDnD/vS2ptl+fkz21kQFu10K/aX/2TH8owP2Z73/feGoLPdNXSWriM72uZx0USmpVWt/lWH9G9ha5JWmo3/wjvfefpHNDvL+962tyKfOF/ThT+M+XVZ713iLAMs0C3Q+h3Ea/5WZ6bu822W9OUiQ1KGH+sNe3Sn/eFxu7pERJbQPM21DSF7LfyPs/1vMU5PVKUmNJe2S/kejqV+vuWf+f+kz7jQK8JnpqCYHux+r0E/UBlPsGlC0E/+GpdwtSXy1pe5AHVa8yjPFfnheI+DDGeG6gFyCdDINuftOf9Ey/zGfaUM+0h4P04X2iXhLOupb94jhZ0jRJf5X9Ql0k6S6/+b0hdWkJ98Nx3yeR7HeqOyUZz//DPW0M9Mw73WfezzxPWN/1ulr21mijAP3FetpeEeB2LfT5P0H2iTJ7Az3BJc0Lct8s9Ez/RYBlMj21X5fUd3nWn6TBnvmzg9z2jaqYEIz04yzQ82WVpEOSGoew/EKVLQQTZW91F/g81otkB91jkjr4zR/2+lYZQrCUMd/raS/Fb3qegrxe+dwfI0p53HX1/O8Nwemhjqs6/dTE3aGhOEf2A/1qY8zVAeoJkpobY5palrXLZ/pR2S+8ARljLpV0u+yzJJup+O7mZrLfEZfHi5IulL3rc4yn3wTZu0i3S/Ld5eQ9DnRakONX3mNvXfyWK02G3/+WpJsty3rBb7q3/9RAx4Vkbw3Eyt7iWeWZ9pGkm2TvUlot6XxJWy3LWmaMWSXPrlFjTHPZ71rftyyrwDOtjuwtiZ2SRpnAV2zky769JekkqbbsLe4DAeqLZO+ODGZlgGmbPb8bl9K3r3DXXx/P9Gz/GS3LKjTGLJK9SyvSyvQ4C/P58qqkRyStM8a8Ifs2LrYsa0ckboAkWZaVL+lWY8wkSRdJGiB7nfaTHRy3GmN+a1nWvz2LVNr6NsZ0kzRa9tZpa0m1/GY5JcBiwV6vvPdXryD311me310krZN9+36SNM6z2/td2btH11iWVRjGzaiS3BqCTWXfdv8Xc3/1ZB8n8Npued4a+TPGjJT9bnGP7P3mm2RvTViyj4v0kv1Os7z+IXvf/h+MMeM9D8Jfy97v/5jlPNbV1PM7UND7qhfOACzLMpJkjKkr+wn1Z0nPGGN+sCzL9ziIt//RYfT/oewQvMAYs1b2rut3fWpjjDENZYejkfN4YGPPtOYq/b4tSUPP75+D1INNlyRZlrU3wGTv/RIbxjjCXX+ljXtbGH2HI+zHWbjPF8uyHjXG7JS9Z+BuSaMkWcaYbEmjLcsK9MajTCzL+ln28feXPGNtIukh2W98njfGtLUs65gqaX0bYwbKfnPoPf79L9mvAUWy3ywOVeDXlmCvV97765ZSuq4nSZZl7feMIVPSZZJ+5anvNMbMkTTV+0a0OnJrCO6TFGNZVpMwlwsWgHGydx1tk9THsqytfvWInZlnWdYRY8zfZD8hL5S0QPZWoeR50vrY5/k91LKsf0VqDD5jOSTpA2PMbyR9KuklY0wny768w7f/hlbxk4yC8YboL3TyoP6HPrXxsoPxAr/5fftbbVlWH5Wdd6wtg9SDTY+0cNefd/5g42tV/iGV2G9Ij7OyPl8sy3pZ0svGPnt7kOzrU2+S9J4xpnMktwr9+t1tjLlN0i8lnSp7D8SnKtv6LvL8Dvba2yjAtDTZeybOsyxroW/BGDNedggGHHqQ6d5x97IsK+ieLUdDlvWjpJuNvXulq+w3oSMkpcu+ymBSKO1URTXuEgkP7yZ6sHfdyyQ19uxiiIRmsh+8SwI8oevp5G6TcMZYkhc9v2/w7Ba8WNJnlmWt8Ztvmed3chn6CJnniTRP9ok795Snf8uytsneBZMse5eUdDIEF8venXmB7CfhHtm7TL3LHpT0paRunnfvZbVe0hFJPY0x9QPUh5SjbX9FKvlxKoW+/j71/E71LxhjYhXZcfsKd5xleb6cYFnWXsuy3rUs6xbZz4UmsncTehV62irLcytYn0Wyj0lK9t4GqWzre4/nd7sAyzTQyV2Rvs6QtNs/AIP1HYIyvy5Yti8ty5ot+024ZG+5V1s1NQT3yH4XdGqQ+izP73kmwLVWxpi6ns3/UG2XvSunr+dJ7G0nXvapzYGu6yttjEFZ9nV538h+B3i77NOcXwww6z9lf7rFCGPMJYHaMvb1THXCHUMAU2UH1P3GGO9xrydlH3udZYwp9uT2XP8U6In4kaQ6so/DfGNZ1mbJ3gqWfdnCb2Ufa1noeXHy9ajsY7rPm8DXezb2HNcIyrOr66+yd3el+S3fS9L1JS0fpl0K8ILoEe76WyJpg6QUY4z/1sGdqpjjgVL4j7Owny/GmPNM4IO8LTy/fT/xx3sII6znljEmwxjTPkjtKkmdZT9vv/BMDnt9e44xr5c02PhcD+kJzUdlb/H5y5PUxBjT029MN+vkrslwvCD7pK8MY0x//6IxJsY4r23sZowJtLXrnVatP22pRu4OtSzroDFmuaRkY8yrkr6W/e7wX5ZlfWZZ1ofGmHGSHpD0jTHmXdlnJdaTfe1LquyTHy4K3EOx/oqMMU/Ivu7pc2PMP2W/EJ8n+13qx56/Qx5jCN2+LPu6t0myjze9GmBcBcaY/5N93dZ/jDFLZF+fdlj2C2+S7MsSWqucD2TLsn4yxjwjO7jGSBpvWdZ6z3Vuz0v60hizwHM742W/QCXLvuaps19zH8p+EWkh+1IH/9q5Pn/7j+N5Y0xf2ceOvjXGvCf7eFMTSafL3mJ4Qfabh5KMk721OcYYM0D2C15r2QH8rux3v/4BXBYfSrrWGDNf9pZFgaQcy7Jywl1/lmVZnhfG/0l6yxjje93aBbJ3nYf0mPZzebBwkH1i0mvhPM7K8nyRfSz8oDFmmexQMJ7bnyT7pKAPfOb9UPbxybc9z+0jkn6wLOuVUm7nPZImG2NWyz65aYfsN0J9ZB/7Pi7pds8JNOVZ3w/LPo6+2Bjzd9knsJwn+35dK/t4qK/HZIfdIs+hkH2yT9YZIvsavatKuV0OlmXt8oT6PyQtM/aHUHwp+015O89tbaqTJ99cKOlhY8xS2Y+/7bL3+gyV/Rx4OJz+q5xon55a3h8Fv3btDNkXde6SfUdZkm70m2eIpL9J2iL7NP8dsp+8j0rq5zdvnqS8EsYRJ/t05XWyn3TbJL0iO1Rf9PTfPtQxKsglEj7Lnio7NC1J80tZRy0kPSj7Hexh2ZcVfCP7CfQHSXHlWdc+9ZaydxkdktTSZ3oPzzr4QfbW4m7PWOZKOj9AO418btvVfrVzvOOQ1KWEsfxaJz8E4Zjn/lghe4u1c4DbtTBAG6fIPs66w3OfrpF9/PUqzzKj/OZfGGz96OT1i/6PwRayL67/2ec2T/abJ9z111f2C/ABz88HnvU2WWW7RKKkn8fK8jhTmM8X2W9a/iH7ouzDnnWwWvYbrvp+446V/aES3+nk5Q7F7t8At3eI7Mt+Fsl+45Qv+7G8Qfbu/h5Blgt7fUu6WXbw5Htu+1zZwRPwMST78bzM0/5e2ZdKpZTwuMpTCa9Xnnnay97b8I3sIN4veyv1FUmX+8zXRfZroveNQb6n/TclDQrlsVSVf7zXYgEIkTFmmuwPU77Isqz3oj0eVE2eyw8yFOCEFlQdNfWYIFBuQY4X95B9iv5uBbg+DED1UiOPCQIRstIYs1H27r1Dsi/6vlT2m8fbrJOfYQqgmiIEgeDmyj4B5jpJ9WUfi3lP0kx2bwE1A8cEAQCuxTFBAIBrlbg71Ph9uSMAANWN5fm840DYEgQAuBYhCABwLUIQAOBahCAAwLUIQQCAaxGCAADXKvcnxgwcOFDDhw9X06ZNFRNDplakoqIi7dq1S3PmzNGyZctKXwAAUKISPzGmtOsEBw4cqPT0dDVv3jziA0NwO3bsUFZWFkEIACGosOsEhw8fTgBGQfPmzTV8+PBoDwMAqr1yhWDTpk0jNQ6EiXUPAOVXrhDkGGD0sO4BoPx4JQUAuBYhCABwLdeH4OTJk3X99deHtcxtt92msWPHnvh/2bJleu211yLS9pYtW5SUlKTc3NywlgMAhM/13yw/bNgwHT16NKxlxo4dq7i4k6tu2bJl+uijj/S73/2u3G0DACqP60Owbdu2YS/ToUOHCmsbAFB52B3qs8ty/vz5SkpK0saNGzVixAglJyfrqquu0kcffeRYxnd36LPPPqtXX31VW7duVVJSkpKSkjR58uRibUvSzp07lZWVpaFDh2rIkCG68sor9fTTT6ugoKBybiwAwMH1W4KBpKWl6YorrtAf//hH/fWvf9XEiRP1zjvvqGXLlsXmHTp0qDZt2qSVK1fq4YcfliQ1btw4YLt79+5VgwYNdM8996h+/fratGmT5s2bpz179mjChAkVepsAAMVVSAi2efZZtZk3r0zLbrnlFm259dYIjyg8v/vd73TZZZdJkjp37qyLLrpIixYt0pVXXlls3pYtW6pZs2ZKSEhQjx49Smz3jDPO0KhRo07836tXL9WuXVtTpkzR6NGjFR8fH9HbAQAoGVuCAQwYMODE340aNVLjxo21ffv2crdrWZZef/11vfPOO9qyZYvy8/NP1LZt26Z27dqVuw8AQOgIwQDq16/v+D8+Pt4RWGX1+uuv64knntD111+vPn36qH79+lq3bp1mzJgRkfYBAOGpkBDccuutUd+lWRV98MEHOv/88x0ffv39999HcUQA4G6uPzs0EkLdUszPz1dCQoJj2oIFCypqWACAUrA7NALat2+v3bt3a/78+erYsaMaNWqkNm3aFJtvwIAB+utf/6pu3bqpbdu2WrBggTZv3hyFEQMAJEIwIn7xi19o5cqVmj17tvbs2aNLL730xLWCvoYNG6Y9e/bomWeekSSdd955uv/++3XvvfdW8ogBAFI5v1n+v//9r5o1axbxQaF0O3fu1MUXXxztYQBAlVdh3ywPAEB1RggCAFyLEAQAuBYhCABwLUIQAOBahCAAwLUIQQCAaxGCAADXIgQBAK5FCEqaN2+eLrnkEvXv31+XXXaZkpKStHHjRklSQUGBnn32WW3YsCHKowQARJrrPzt03bp1evbZZzV8+HD17dtXtWrVUn5+vtq2bSvJDsF58+apdevW6tSpU5RHCwCIJNeHYF5eniTp6quvVr169aI7GABApXL17tDJkycrIyNDkv2NDklJSVq1apVjd2hqaqokKSsrS0lJSUpKStKWLVuiNmYANZNlWRXyg5K5ektw2LBhatmypZ5//nk9/fTTSkxM1O7dux3zPP3007rjjjt00003aciQIZLEN2cAQA3h6hBs27btiWN/Xbt2VZ06dbRq1SrHPF27dj0xb48ePSp9jACAihPxEEz6T1KkmwzbJ5d+Eu0hAACqAVdvCQJAtFTW8Tr/fowJ+v2yruTqE2MAAO4W8S3BmrYrMj4+XpJ07NixKI8EABBp7A4tRXx8vNq0aaMPPvhAHTt2VEJCgs4888wT4QgAqL7YHRqC8ePHa+/evRoxYoRuuOEG7dixI9pDAgBEgCnp4KwxpsQjt//973+5Zi5Kdu7cqYsvvjjawwBQRtG6kN2NJ8ZYlhX0RrMlCABwLUIQAOBahCAAwLUIQQCAaxGCAADXKlcIFhUVRWocCBPrHgDKr1whuGvXrkiNA2Fi3QNA+ZUrBOfMmcOF41GwY8cOzZkzJ9rDAIBqr1wXy0vSwIEDNXz4cDVt2lQxMRxirEhFRUXatWuX5syZo2XLlkV7OADKgYvlK09JF8uXOwQBAOEjBCsPnxgDAEAAhCAAwLUIQQCAaxGCAADXIgQBAK5FCAIAXCsu2gMAADfyv1Shoi6ZcOMlEeFgSxAA4FqEIADAtQhBAIBrEYIAANfixBgAqAI4gSU62BIEALgWIQgAcC1CEADgWoQgAMC1CEEAgGsRggAA1yIEAQCuRQgCAFyLEAQAuBYhCABwLUIQAOBahCAAwLUIQQCAaxGCAADXIgQBAK5FCAIAXIsQBAC4FiEIAHAtQhAA4FqEIADAtQhBAIBrEYIAANciBAEArkUIAgBcixAEALgWIQgAcC1CEADgWoQgAMC1CEEAgGsRggAA1yIEAQCuRQgCAFyLEAQAuBYhCABwLUIQAOBahCAAwLUIQQCAaxGCAADXIgQBAK5FCAIAXIsQBAC4FiEIAHAtQhAA4FqEIADAtQhBAIBrEYIAANeKC2dmy7KKTTPGRGwwAABUJrYEAQCuRQgCAFyLEAQAuBYhCABwrRJPjAl0Ikxp83CiDACgumBLEADgWoQgAMC1CEEAgGsRggAA1yIEAQCuRQgCAFyr1M8Onb18tg4XHNaI/iNUL6FeZYwJAIBKYUq6FnB//n6r/WPttefoHjWt3VT3nXOfRvQfoQaJDULvgOsGAQBRZFlW0CAqcXfo86uf156jeyRJu47s0oSPJqj9Y+01NWeq9h3dF+FhAgBQuUrcEjxWeMx6ee3LmpY7TXl78xy1RrUaadSAURo5cKQa1WoUvAO2BAEAUVTSlmCJISjJkqSCwgK98tkrmpY7Td/t+c4xQ4PEBho5YKRGDRylJrWbFO+AEAQARFG5Q9CroLBAr37+qqblTtPG3RsdM9ZPqK+7B9ytewbeo6Z1mp7sgBAEAERRmUPQGBO4GCOpu6QUSc2cpXoJ9XRn0p26b9B9alanWYCFA/ZT6jx8UDcAoCwiH4InZpBeW/uapuRM0Vc7v3KU6sbX1YikEbpv0H1qUbdFyc0QggCAClJxIWg3rsKiQr257k1NyZmiL3d86ajXia+jO/rdofsH3a9W9VoF66e0bghBAECZVHgIehVZRXr7q7eVlZ2lz7d/7pivVlwt3d73do0ZPEat67f276e0bghBAECZVGgIBuhMRVaR3ln/jrKys7T257WOemJsom7te6vGDh6rUxqcEm7zYSEoAQCVHoJeRVaR5m+Yr8zsTK3ettoxX0Jsgob1HqZxQ8apXcN24XYTEkIQABC1EPSd9u+v/62snCyt3LLSUUuITdBNZ9+k8cnjdWrDU8PtrkSEIAAg6iHoW1uwcYEyszO1/Kfljlp8TLxuPPtGTUieoPaN2ofbbUCEIACgyoSg7zzvf/u+MrMztfTHpY5aXEycbuh1gyYkT1CHxh3C7d6BEAQAVGoIBui8xNqH33+ozOxMLdq0yFkskrRWUq6k3aGFa2kIRQBwnyobgr7zfJz3sTKzM5XzQ46zWCTpM2nDcxt0VtOzyjUWQhAA3KfKh6Cv7LxsnTv5XOl05/QYE6Prul+ntJQ0dW7WuUxjIQQBwH2qVQh6+pVOlZQqqaNfTUbXdL9Gk1ImqWvzruG3CwBwlaiGYJABhTzvks1LlJmdqfe/fd8x3cjo6m5XKy05TT1a9ojY2AhKAKhZqnUIepl2xt4yPLN47couVyo9NV09W/Ys99gIQQCoWWpGCHrDqY3sMOxUfJ7LO1+u9JR09W7du8xjIwQBoGYpKQRjKnMgEbFF0uvSqltXaWinoY7SO+vfUZ9n+2joG0O1asuq6IwPAFBtRGVL0F+5tgxbyd4y7FJ8nkvPvFTpqenqf0r/Mo+NLUMAqN6q3O5Qf+UKQa+Wsr/pvlvxeS864yJlpGZoYNuB5e8HAFCtuCMEvVroZBj6zXJhhwuVkZqhwacOLn8/AIBqwV0h6NVcuu7p6/TGF2/IkrP9C06/QBmpGUo+Lbn8/QAAqjR3hqCn3fU712ta7jS99vlrKrKKHPXz2p+n9NR0ndv+3HL1AwCouqp8CPoLNKZyh1FTScmSeqrYObEpp6UoPSVd559+fkj9EIwAUH0Qgr4ayz5m2EvFwnBwu8FKT03XhR0uLLE/QhAAqg9CMJDGkoZIcUlxOl503FEa2HagMlIz9KuOvwrYLyEIANUHIViC7/d8rwdyH9ALa15QQVGBo9b/lP5KT0nXJWde4uifEASA6qPahWBFKTFcG0oaIiUMTNCxwmOOefq16af0lHT9+qxfs2UIANUMIegRyhbm5n2b9dCihzTv03nKL8x31Hq36q301HQN7TSULUMAqCYIQY9QQtA7z5YDWzRj8QzNXTVXR48fdczTq2UvTUqZpCu6XKEYE0MIAkAVRgh6hBOCXlsPbNXDSx7WMyuf0ZHjRxy1Hi16aFLKJP22x2+lGrWmAKDmcG0IRuIifG8bPx/8WTOXzNSclXN0uOCwY56uzbtqUsokXd31asXGxJapHwBAxSAEw1DaluGOQzv0yNJH9OSKJ3Wo4JCj1rlZZ6Ulp+na7teWGoaEIABUDkIwDKWFoNfOwzs1a+kszV4xWweOHXDUzmp6ltKS03Rdj+sUFxMXUj8AgIpBCIYh1BD02n1kt2YtnaUnVjyh/fn7HbUzmpyhickT9YeefygWhoQgAFQOQjAM4Yag154je/T48sf12LLHtC9/n6PWoXEHTRgyQdf3ul7xsfEB+wEAVAzXhqC/SH0STUnrbO/RvXpi+ROatWyW9h7d66i1b9Re44eM141n36iE2ISw+w0loAlXAHAiBD0qIwS99h3dpydXPKlHlz2q3Ud2O2qnNjxV44eM15/O/pMS4xJD7pcQBIDwEYIelRmCXgfyD+ipT57SzCUztevILketbYO2Gjd4nG7uc7NqxdUqtS1CEADCRwh6RCMEvQ4eO6g5n8zRzCUztePwDketTf02Gjd4nIb1Gaba8bWDtkEIAkD4CMEKFk4oHjp2SM+sfEYzlszQ9kPbHbXW9VprzOAxurXvraoTXyciYyMUAbgdIVjByrJleLjgsOaunKsZS2Zo28FtjlrLui01etBo3d7vdtVNqFuusRGCANyOEKxgZQlBryMFR/Tcp8/pwcUPasuBLY5a8zrNNXrQaN2RdIfqJdQrU/uEIAC3IwQrWHlC0Ovo8aN6fvXzemDRA/px/4+OWrM6zXTfOfdpRNII1U+sH1a7hCAAtyMEq6Bg6z3/eL5eWPOCHlj0gDbt2+SoNandRPcOvFd3DbhLDRIbhNQPIQjA7QjBKqi0rcdjhcf00pqXNH3RdOXtzXPUGtdqrHsG3qO7BtylRrUaldgOIQjA7QjBKijUXagFhQV6ee3LmpY7Td/v/d5Ra5jYUKMGjtLIASPVuHbjgMsTggDcjhCsgsI9jlhQWKBXP39VU3Om6ts93zpqDRIb6O7+d+uec+5Rk9pNHDVCEIDbEYJVUFlPpjledFyvf/66puZO1de7vnbU6iXU013979K959yrZnWaSSIEAYAQrKZKum8Kiwr1xhdvaGruVK3fud5RqxtfV3f2v1P3nXOfmtdtHna/BCeAmoQQrKZC2VosLCrU39f9XVNypmjdjnWOWp34Ohreb7juH3S/WtZrGXK/hCCAmoQQrKbC2WVaZBXprXVvKSsnS19s/8JRqx1XW7f3u11jBo9Rq3qtSm2LEARQkxCC1VRZjhsWWUX6x1f/UFZOlj77+TNHrVZcLd3W9zaNGTxGbeq3CdoGIQigJiEEawj/+ypQWHnnKbKK9M/1/1RWTpbWbFvjmCcxNlG39LlFY4eMVdsGbcs0Fr7RAkB1QQjWEOGEoO//87+er6zsLK3auspRS4hN0LDewzRuyDi1a9gurLEQggCqC0KwhihLCPpOf/ebd5WZnalPtnziqMXHxOum3jdp/JDxOq3RaSGNhRAEUF0QgjVEeULQt/7et+8pMztTy35c5qjFxcTpxl43akLyBJ3e+PQS2yEEAVQXhKCLhHoyjWVZ+uC7D5SZnanFmxc7anExcbq+5/WakDxBHZt0jNjYCEUA0UAIuki4Z5RalqWP8z5WZnamcn7IcdRiTaz+0PMPmpg8UWc2PbPcYyMEAUQDIegi5fluw4V5C5WZnamFeQsd02NMjH7X43dKS05Tp2adytw+IQggGghBF4nEF/zm/JCjrOwsffj9h47pMSZG13a/VmnJaerSvEvY7RKCAKKBEHSRSISg1+JNi5WZnan/ffc/x3Qjo992+60mpUxStxbdQm6PEAQQDYSgi4VyRmlpyyzdvFRZOVlasHFBsXmv6nqVJqVMUs+WPYvVSuuLM0oBVAZC0MUiEYJeK35aoazsLP3nm/8Uq13R+Qqlp6br7FZnh9wXIQigMhCCLhbJEPRauWWlsrKzNP/r+cVql3W6TOkp6erbpi8hCKBKIARdrCJC0Gv11tXKysnSO+vfKVb79Vm/1r/v/7e0Jbx+CEEAkUYIIuJ8Hzdrt63V2SPPlroWn+/iMy5WRmqGBrQdUKZ+CEUA5UUIIuICbmG2kJQqOwz9HnK/6vgrpaema1C7QWH1QwgCKC9CEBFX4m7W5pJSJNPDyJJzvl90+IUyUjM05NQhIfVDCAIoL0IQERfKscZ129dpau5UvfHFGyqyihy189qfp4zUDKW2Ty2xH0IQQHkRgqhwJT2ONuzcoOmLpusvn/2lWBimnpaqjNQMndv+3JACryyhyAk4gLsRgqhwoZxRunH3Rk3LnaZX1r6iQqvQUUs+NVnpqem64PQLSgwoQhBAuAhBVLhwPq7tuz3faXrudL209iUdLzruqA1qN0gZqRm6sMOFAYOKEAQQLkIQFa4sn1matzdPD+Q+oBfWvKCCogJHbcApA5SRmqGLzrjIEViEIIBwEYKocOX54O4f9v6ghxY/pD+v/rOOFR5z1JLaJCk9NV2XnnmpjDGEIICwEYKIinCDcfO+zXpo8UOa9+m8YmHYt3VfTUqZpMu7XF7ufgIhFIGaixBEVJQ1nH7a/5NmLJ6huavmKr8w31ncKilb0gbJewkiIQigJIQgoqK84bT1wFbNWDxDz6x6RkePH3UWt0nKkfSVZBURggCCIwQRFZH6gt9tB7dp5pKZmvPJHB05fsRZ3C69cccbuqrrVYqNiS1zH4QgUHMRgqgSQvmUmZIej9sPbdfMJTP11CdP6XDBYUetS7MuSktJ0zXdrgkpDAk9wD0IQVQJ5Q1Brx2HdujRpY/qyU+e1MFjBx21Tk07KS0lTdd2v1ZxMXFB2yAEAfcgBFElRCoEvXYd3qVHlz6q2Stm68CxA47amU3O1MTkifp9z98HDENCEHAPQhBVQqRD0Gv3kd16bNljenz549qfv99R69i4oyYmT9Qfev5B8bHxJfYNoGYiBFFthROKe4/u1ePLHtdjyx/T3qN7ncU9knIlrZVUGFogA6gZCEFUW2XZMtx3dJ9mr5itR5c+qj1H9ziLeyXlSvnL85UQm3BiMiEI1FwlhWBMKA1kyL4uuSw/GeUYOP3Sr4wJ/DN5ctB2G9ZqqLSUNOWNytP086dLvieSNpL0G+mM0Yl6OskoP85ur6rcXvqlX/qNXL+h9BNSCALVUYPEBhqfPF56TNL/JB06WdvcUBr+a6njSOnJ/tLR4CeSAqjBCEHUfMckLZYdhu9LdX2uqvipgXTXJVLHu6UnBkhHCEPAVTgmiGqlLN8I4b/M4YLDmrtyrmYsmaFtB7c5aq3qtdKYQWN0W7/bVCe+TontchwRqB44MQY1RiRC0OtIwRE9u+pZPbT4IW09uNVRa1G3hUYPGq07+t2hugl1Ay5PCALVAyGIGiOSIeh19PhRPffpc3pw0YP66cBPjlqzOs10/zn3a3jScNVPrB9WvwCqBkIQNUZFhKBX/vF8Pb/6eT2w6AFt3r/ZUWtau6nuPede3dn/TjVIbBBSvwCqBkIQ8FPS4z7/eL5eXPOipi+ark37NjlqjWs11r3n3Ku7+t+lhrUaOmqEIlA1EYKAn1C2Do8VHtNLa17S9EXTlbc3z1FrVKuRRg0YpZEDR6pRrUaSCEGgqiIEAT/hfBJNQWGBXvnsFU3Lnabv9nznqDVIbKCRA0Zq1MBRalqnaaSHCSACCEHAT1k+ju140XG9+tmrmpo7VRt3b3TU6ifU14EPDkhLJR0JvDyA6CAEAT/l+db740XH9cYXb2hqzlRt2LXBWcyXtEJ2GB4OsDCASkcIAqUoSygWFhXqb1/+TVNypuirnV85anXj62p40nDdP+h+tajb4sR0jhsClY8QBEpRni3DwqJCvbnuTU3JmaIvd3zpqNWJr6M7+t2h0YNGq2W9loQgEAWEIFCK8oSgV5FVpLe/eltZ2Vn6fPvnjlrtuNq6re9teuzqx6SDgZcHUDEIQaAUkQhBryKrSP9c/09l5WRpzbY1zmKBpFWyP9D7QMS6BFACQhCIgHCD0rIszf96vjKzM/Xp1k8dtcTYRA3rM0zjhoxT2wZtyz02drMCwRGCQASUdWvRsiz955v/KDM7Uyu3rHTUEmITdNPZN2l88nid2vDUMo+NEASCIwSBCCjvLlPLsrRg4wJlZmdq+U/LHbX4mHjdePaNGj9kvE5vfHrYbROCQHCEIBABkTpuaFmW3v/2fWVmZ2rpj0sdtbiYOF3f83pNTJmoDo07hNwmIQgERwgCFSASW4Yffv+hMrMztWjTIkct1sTqj73+qInJE3VGkzPCbptQBE4iBIEKEMktw4/zPlZmdqZyfshx1GJMjH7f4/dKS0nTWU3PCrlNQhA4iRAEKkAkL6vwys7LVmZ2pj7O+9gxPcbE6Lru1yktJU2dm3UutR1CEDiJEAQqQEWEoFfuD7nKysnSB9994JhuZHRN92s0KWWSujbvGnR5QhA4iRAEKkBFhqDXks1LlJWdpfe+fc8x3cjo6m5Xa1LKJHVv0b3YcoQgcBIhCFQRZQ3O5T8uV1ZOlt795t1itSu7XKm37n5L+rnkfsoSjP7tEK6ojghBoIoo79bjyi0rlZWdpflfzy9e/EpStqRthCDgixAEqohI7UL9dOunysrO0j83/LN4cYO0cuZK9W3T1zGZEIRbEYJAFRHp44hrtq3RlJwpevurt4vVLj3zUqWnpqv/Kf0lEYJwL0IQqOHWblurKTlT9Na6t2TJ+bS96IyLlJGaoYFtBzqm+wdaWQKaUER1QAgCNZz3efzl9i81JWeK/vbl34qF4YUdLlRGaoYGnzpYEiEI9yAEgRrO/3n81Y6vNDV3qt744g0VWUWO2vmnn6+M1Ayltk8tsY1QEIKoDghBoIYL9jzesHODpuZO1Wufv1YsDJUnaaHndwltlIQQRHVACAI1TLiB9c2ubzQtd5r+8tlfVGgVOmrJpyYrIzVD559+fqmhRuihOiIEgRqmrGeZfrv7W03Pna6X1r5ULAwHtxus9NR0XdjhwqBhRwiiOiIEgRqmvJdafL/ne03Pna4X176o40XHHbWBbQcqPSVdF51xUbHQIwRRHRGCQA0TqesN8/bm6cFFD+r51c+roKjAUet/Sn+lp6TrkjMvORF+hCCqI0IQqGEifdH9pn2b9OCiB/Xn1X/WscJjjlrf1n2Vnpqu35z1G8XExES0X6AyEIJADRfKJ7uEEpw/7v9RDy16SPM+naf8wnxHrXer3kpPTdfQTkNL3DLkU2ZQ1RCCQA0XqRD02nJgi2YsnqG5q+bq6PGjjlqvlr00KWWSruhyhWJjYss0FqAyEYJADRfpEPTadnCbHl78sJ5e+bSOHD/iqHVv0V1fPPWFtE7y/XAaQhBVDSEI1HAVFYJePx/8WTOXzNSclXN0uOCws7hdUo6kLyVZhCCqHkIQQEROptlxaIceWfqInlzxpA4VHHLUOjfrrLTkNF3b/dqAu0lLUtrnmBKkKA9CEEBEzyjdeXinZi2dpdkrZuvAsQOO2llNz1Jacpqu63Gd4mLiQmqPEERFIgQBRPyyCknafWS3Zi2dpSdWPKH9+fsdtTOanKGJyRP1+x6/V3xsfIntEIKoSIQggAoJQa89R/bo8eWP67Flj2lf/j5HrUPjDpowZIKu73V90DAkBFGRCEEAxZQlaEoL0r1H92r28tmatWyW9hzd46id1vA0TUieoBvPvlEJsQnhD9gPwYlQEYIAiqmIEPTan79fs5fP1qPLHtXuI7sdtXYN2mn8kPG6qfdNSoxLDH3AfghBhIoQBFBMRYag14H8A3rqk6c0c8lM7Tqyy1Fr26Ctxg0ep5v73KxacbXCalciBBE6QhBAMZURgl4Hjx3UnE/maOaSmdpxeIej1qZ+G40dPFa39LlFteNrh9wmIYhQEYIAiqnMEPQ6dOyQnln5jGYsmaHth7Y7aq3qtdLYwWN1a99bVSe+TqltEYIIFSEIoEIEev0oLZwk6XDBYc1dOVczlszQtoPbHLWWdVtq9KDRur3f7aqbUDei4yUY3YkQBFAhyhqCXkcKjui5T5/Tg4sf1JYDWxy15nWaa/Sg0boj6Q7VS6gXkfESgu5ECAKoEOUNQa+jx4/q+dXP64FFD+jH/T86as3qNNN959ynEUkjVD+xfrnGSwi6EyEIoEJEKgS98o/n68U1L2r6ounatG+To9akdhPdO/Be3TXgLjVIbFCm8RKC7kQIAqg0kfhkmmOFx/TSmpc0fdF05e3Nc9Qa1Wqkewbeo7sH3K1GtRoFbaOs36RBUNY8hCCAShPJj2crKCzQK5+9omm50/Tdnu8ctYaJDTVq4CiNHDBSjWs3LrYsIQgvQhBApamIzygtKCzQq5+/qqk5U/Xtnm8dtQaJDXR3/7t1zzn3qEntJiemE4LwIgQBVJqK/KDu40XH9frnr2tq7lR9vetrR61eQj3d1f8u3XvOvWpWpxkhiBMIQQBRE4mL2P3bKCwq1BtfvKGpuVO1fud658zHJK2Qtr+zXc3rNg+7L3+EYvVHCAKImooIQa/CokL9fd3fNSVnitbtWOeo1Ymvo+H9huv+QferZb2WYffpRQhWf4QggKipyBD0KrKKFNsjVkqR5Jd3teNq6/Z+t2vM4DFqVa9V2H0TgtUfIQggaiojBE+0ayR1lpQqyS/vasXV0m19b9OYwWPUpn6bkPsmBKs/QhBA1FRqCJ74R3pn3TvKzM7U6m2rHfMlxibqlj63aOyQsWrboG147aJaIgQBuI5lWbIsS//++t/KysnSyi0rHfWE2ATd3PtmjRsyTqc2PDXkdkM565TgrFoIQQCu4/vaZlmW/rvxv8rMztSKn1Y45ouPiddNvW/S+CHjdVqj00ptlxCsfghBAK4T6LXNsiy99+17ysrO0tIflzpqcTFxurHXjZqQPEGnNz49aLuEYPVDCAJwnZJe2yzL0gfffaDM7Ewt3rzYUYs1sbq+1/WamDxRHZt0LLYsIVj9EIIAaryyfFKNZVn6OO9jZWZnKueHHEct1sTqDz3/oInJE3Vm0zPLNTaCM7oIQQA1Xnk/rm1h3kJlZWfp47yPHdNjTIx+1+N3SktOU6dmncrUNiEYXYQggBovUp9ZmvtDrjKzM/Xh9x86pseYGF3T7RpNSpmkLs27hNUmIRhdhCCAGi/SH9y9eNNiZeVk6f1v33dMNzL6bbffKi0lTd1bdA+pLUIwughBAK4Tyrfeh7Lc0s1LlZWTpQUbFxSb98ouVyo9NV09W/YssZ+yBDShGDmEIADXiVQIeq34aYWysrP0n2/+U6x2RecrlJ6arrNbnR2wH0IwughBAK4T6RD0WrVllbJysvSvDf8qVrus02VKT0lXv1P6hdVmIIRg5BCCAFynokLQa/XW1crKydI7698pXtwgKVvSlvDa9EUIRg4hCMB1KjoEvdZuW6spOVP01ldvFS9+IylbsjYTgtFECAJAhAQN1xayv8+wm+yvdPLxy46/VEZqhga1GxS0XUKv4hCCABAhpW5hNpeue/o6vfHFG7LknPcXHX6hjNQMDTl1SMltIKIIQQCIkFB2s1qWpfU712ta7jS99vlrKrKKHPXz2p+njNQMpbZPDdoGIocQBIAICTUEvb7e9bWm5U7Tq5+9qkKr0DFf6mmpSk9N13ntz1NMTEzFDBiEIABUpFBOptm4e6Om507Xy2tfLhaGQ04doozUDF1w+gVhbRGWZevRjZ9UQwgCQAUK54zS7/Z8pwdyH9CLa1/U8aLjjtqgdoOUnpKuX3b8ZUjhRAiGhhAEgApUlusA8/bm6YHcB/TCmhdUUFTgqA04ZYDSU9N18RkXlxhShGBoCEEAqEDl+fDuTfs26cFFD+rPq/+sY4XHHLV+bfopPSVdvz7r1wHDihAMDSEIAJWoLKG4ed9mPbT4Ic37dF6xMOzdqrfSU9M1tNPQSgmtmhaMhCAAVKLybBn+tP8nzVg8Q3NXzVV+Yb6j1qtlL6WnpuvyzpcrxlTc2aSEoLdICAJA2CLx3YZbD2w9EYZHjh9x1Hq06KFJKZN0ZdcrKyQMCUFvkRAEgLBF8gt+tx3cpplLZmrOJ3OKhWG35t00KWWSrup6lWJjYiPWJyHoLRKCABC2SH/LvSRtP7Rdjyx5RE998pQOFRxy1Lo066K0lDRd0+2aiIQhIegtEoIAUClCDc6dh3fqkSWP6MlPntTBYwcdtU5NOyktJU3Xdr9WcTFxERtbdQ9FQhAAqrhwtx53Hd6lWctm6YnlT+jAsQOO2plNztTE5In6fc/fRyQMCUEAQIUq6y7U3Ud26/Flj+vx5Y9rX/4+R61D4w6amDxRf+z5R8XHxpd5bIQgAKBClfc44t6je/XE8ic0a9ks7T2611Fr36i9JgyZoBvOvkEJsQlht00IAgCqnECv3/uO7tPsFbP16NJHtefoHkfttIanafyQ8fpT7z+VGIalhV4o36RRlRCCAFADlfT6vT9/v55a8ZQeWfqIdh3Z5ai1a9BO44aM0829b1ZiXGKxZQlBb5EQBIAqK5RdqAePHdScT+bo4SUPa+fhnY7aKfVP0bgh4zSszzDViqt1Yjoh6C0SggBQZYVzHPHQsUN6euXTenjJw9p+aLuj1rpea40dPFa39r1VteNrE4InioQgAFRZZTmZ5nDBYc1dOVczlszQtoPbHLWWdVtqzOAxui/1PqkgSANB+iUEAQBVTtCvToqT1FfSEEn1ncu0qNtC959zv+5IukP1EuqVqd+qFIqEIAC4VKnfHxgnqY/sMGzgLDWr00z3nXOfRiSNUP1Ev6QsBSEIAIi6kL9EN1ZSb0nJkho6S01qN9F959ynO/vfqQaJDQIsXBwhCACIurC/ST5WembZM5q+aLo27dvkKDWu1Vj3DLxHdw+4Ww1rNQzSQIj9VCJCEAAQULAMOFZ4TC+vfVnTcqcpb2+eo9aoViONGjBKIweOVKNajQIuH0oIhh3QZUQIAgACKu0M04LCAv3ls79oWu40fbvnW0etQWIDjRwwUqMGjlKT2k0cNUIQAFDlhXqZxfGi43rt89c0NWeqvtn9jaNWP6G+7up/l+495141rdNUEiEIAKgGwr3W8HjRcf31i79qSs4Ubdi1wVGrl1BPdybdqXvPuVct6rUIu29CEAAQdaEEY2FRof725d80NXeq1u1Y56jVja+r4UnDdf+g+9WibulhWJJIBCMhCAAIWThbh0VWkd5c96aysrP05Y4vHbXacbV1R787NHrwaLWq16pMYyEEAQCVqiwfx1ZkFentr97WlJwp+uznzxy1WnG1dHvf2zVm8Bi1rt86rHYJQQBApSrPF/wWWUX65/p/KisnS2u2rXHUEmMTdWvfWzV28Fid0uCUkNojBAEAlaq833LvbeNfG/6lrJwsfbr1U0ctMTZRw/oM09jBY9WuYbsS2yEEAQBRFcpZnMGyxLIs/eeb/ygzO1Mrt6x01OJj4nVz75s1bsg4ndboNM4OBQBUPeUJQd/6go0LlJmdqeU/LXfU4mLidGOvG/XcDc9Je8s93EB9E4IAgLKJRAj6zve/7/6nzOxMLdm8xFkslLRWUq6kPWUba5A+CUEAQNlEMgR95//w+w+VmZ2pRZsWOYtFkj6TlCNpd3hjDdIXIQgAqDyhhOKJMG0v6VzPbx8xJka/7/F7paWk6aymZ5XcRsl9E4IAgMoTVgh6LPx+oTKzM/Vx3seO6TEmRtd2v1ZpyWnq0rxLiW0E6ZsQBABUnrKEoHeZ3B9yNSVniv733f+c88vomu7XKC05Td1adAvYRpC+CUEAQOUpTwh6Ld28VFk5WVqwcYFzORld1fUqTUqZpJ6teobSNyEIAKg4kbjAPpjlPy5XVk6W3v3m3WK1/+vyf0pPSVevVr1KaoIQBABUnIoMQa+VW1YqKztL87+eX6x2eefLlZ6Srt6tewdalBAEAFScyghBr0+3fqopOVP0zvp3itV+c9ZvlJ6arn5t+vlOJgQBABWnMkPQa822NZqSM0Vvf/V2sdolZ16ijNQM9T+lv0QIAgAqUjRC0Ouznz/T1JypenPdm7LkHMfNvW/Wc5c9FzQEYyp8dACAGs8Y4/gJZZ5QlglFz5Y99ber/6bP7/hc13a/VsZnw69P6z4lj5stQQBApAXKltKCLlJbk1/t+EpTc6cq94dcfX3X16oVV4vdoQCAyhPNEPQ6dOyQ6ibUlTgmCACoTFUhBH27DlaIq6geAQDuVZZjfGVZJpSwLSlcOTEGAOBahCAAwLUIQQCAa3FMEABQbZXn+kKJLUEAgIsRggAA1yIEAQCuRQgCAFyLEAQAuBYhCABwLUIQAOBahCAAwLUIQQCAaxGCAADXIgQBAK5FCAIAXIsQBAC4FiEIAHAtQhAA4FqEIADAtQhBAIBrEYIAANciBAEArkUIAgBcixAEALgWIQgAcC1CEADgWoQgAMC1CEEAgGsRggAA1yIEAQCuRQgCAFyLEAQAuBYhCABwLUIQAOBahCAAwLUIQQCAaxGCAADXIgQBAK5FCAIAXIsQBAC4FiEIAHAtQhAA4FqEIADAtQhBAIBrEYIAANciBAEArkUIAgBcixAEALgWIQgAcC1CEADgWoQgAMC1CEEAgGsRggAA1yIEAQCuRQgCAFyLEAQAuBYhCABwLUIQAOBahCAAwLUIQQCAaxGCAADXIgQBAK5FCAIAXIsQBAC4FiEIAHAtQhAA4FqEIADAtQhBAIBrEYIAANciBAEArkUIAgBcixAEALgWIQgAcC1CEADgWoQgAMC1CEEAgGsRggAA1yIEAQCuRQgCAFyLEAQAuBYhCABwLUIQAOBahCAAwLUIQQCAaxGCAADXMpZlRXsMAABEBVuCAADXIgQBAK5FCAIAXIsQBAC4FiEIAHAtQhAA4Fr/Dz1Ca5n3jQZKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax = plt.subplots(1,1, figsize=(8,8))\n",
    "ax.imshow(im,cmap='gray')\n",
    "l1 = ax.plot(x,y_init,'-r',lw=3, linestyle='--', dashes=(5,5), label='initial')\n",
    "l2 = ax.plot(x,y_irls,'-g',lw=3, label='fit')\n",
    "\n",
    "ax.set_aspect('equal')\n",
    "ax.axis('off')\n",
    "ax.set_title('Iterative Reweighted Least Squares', fontsize=20)\n",
    "ax.legend(loc='upper left', fontsize=15)"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}