bundleAdjustment.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bundle Adjustment using scipy\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We have a set of points in real world defined by their coordinates $(X, Y, Z)$ in some apriori chosen \"world coordinate frame\". We photograph these points by different cameras, which are characterized by their orientation and translation relative to the world coordinate frame and also by focal length and two radial distortion parameters (9 parameters in total). Then we precicely measure 2-D coordinates $(x, y)$ of the points projected by the cameras on images. Our task is to refine 3-D coordinates of original points as well as camera parameters, by minimizing the sum of squares of reprojecting errors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's start solving some real bundle adjusment problem. We'll take a problem from http://grail.cs.washington.edu/projects/bal/."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import urllib\n",
    "import bz2\n",
    "import os\n",
    "import numpy as np\n",
    "import open3d as o3d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First download the data file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "BASE_URL = \"http://grail.cs.washington.edu/projects/bal/data/\"\n",
    "\n",
    "DATASET_NAME = \"ladybug/\"\n",
    "FILE_NAME = \"problem-49-7776-pre.txt.bz2\"\n",
    "# FILE_NAME_BIG = \"problem-412-52215-pre.txt.bz2\"\n",
    "\n",
    "# DATASET_NAME = \"dubrovnik/\"\n",
    "# FILE_NAME = \"problem-88-64298-pre.txt.bz2\"\n",
    "\n",
    "URL = BASE_URL + DATASET_NAME + FILE_NAME"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.isfile(FILE_NAME):\n",
    "    urllib.request.urlretrieve(URL, FILE_NAME)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now read the data from the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_bal_data(file_name):\n",
    "    with bz2.open(file_name, \"rt\") as file:\n",
    "        n_cameras, n_points, n_observations = map(\n",
    "            int, file.readline().split())\n",
    "\n",
    "        camera_indices = np.empty(n_observations, dtype=int)\n",
    "        point_indices = np.empty(n_observations, dtype=int)\n",
    "        points_2d = np.empty((n_observations, 2))\n",
    "\n",
    "        for i in range(n_observations):\n",
    "            camera_index, point_index, x, y = file.readline().split()\n",
    "            camera_indices[i] = int(camera_index)\n",
    "            point_indices[i] = int(point_index)\n",
    "            points_2d[i] = [float(x), float(y)]\n",
    "\n",
    "        camera_params = np.empty(n_cameras * 9)\n",
    "        for i in range(n_cameras * 9):\n",
    "            camera_params[i] = float(file.readline())\n",
    "        camera_params = camera_params.reshape((n_cameras, -1))\n",
    "\n",
    "        points_3d = np.empty(n_points * 3)\n",
    "        for i in range(n_points * 3):\n",
    "            points_3d[i] = float(file.readline())\n",
    "        points_3d = points_3d.reshape((n_points, -1))\n",
    "\n",
    "    return camera_params, points_3d, camera_indices, point_indices, points_2d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "camera_params, points_3d, camera_indices, point_indices, points_2d = read_bal_data(FILE_NAME)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "camera_params: (49, 9);\n",
      "points_3d: (7776, 3);\n",
      "camera_indices: (31843,); \n",
      "point_indices: (31843,); \n",
      "points_2d: (31843, 2)\n"
     ]
    }
   ],
   "source": [
    "print(f\"camera_params: {camera_params.shape};\\npoints_3d: {points_3d.shape};\\n\"\n",
    "        f\"camera_indices: {camera_indices.shape}; \\npoint_indices: {point_indices.shape}; \\n\"\n",
    "        f\"points_2d: {points_2d.shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we have numpy arrays:\n",
    "\n",
    "1. `camera_params` with shape `(n_cameras, 9)` contains initial estimates of parameters for all cameras. First 3 components in each row form a **rotation vector**, next 3 components form a translation vector, then a focal distance and two distortion parameters.\n",
    "2. `points_3d` with shape `(n_points, 3)` contains initial estimates of point coordinates in the world frame.\n",
    "3. `points_2d` with shape `(n_observations, 2)` contains measured 2-D coordinates of points projected on images in all the observations.\n",
    "4. `camera_ind` with shape `(n_observations,)` gives the index of the camera (from 0 to `n_cameras - 1`) associated with a particular observation.   \n",
    "    - TODO-for-students: Basics - How do we know this (`camera_ind`) information in practical setting? Explain.     \n",
    "5. `point_ind` with shape `(n_observations,)` contains indices of 3D points (from 0 to `n_points - 1`) involved in each observation.\n",
    "    - TODO-for-students: Basics - How do we know this (`point_ind`) information in practical setting? Explain.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 0 0 0 0 0 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2] \n",
      " [ 0  1  3 26 29 36  0  1  4  8 20 38 47  0  1  2  3  4  5  6  7  8  9 10\n",
      " 11 12 13 15 17 20]\n"
     ]
    }
   ],
   "source": [
    "print(point_indices[:30],\"\\n\",camera_indices[:30])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n_cameras: 49\n",
      "n_points: 7776\n",
      "Total number of parameters to estimate: 23769\n",
      "Total number of residuals: 63686\n"
     ]
    }
   ],
   "source": [
    "n_cameras = camera_params.shape[0]\n",
    "n_points = points_3d.shape[0]\n",
    "\n",
    "n = 9 * n_cameras + 3 * n_points\n",
    "m = 2 * points_2d.shape[0]\n",
    "\n",
    "print(\"n_cameras: {}\".format(n_cameras))\n",
    "print(\"n_points: {}\".format(n_points))\n",
    "print(\"Total number of parameters to estimate: {}\".format(n))\n",
    "print(\"Total number of residuals: {}\".format(m))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We chose a relatively small problem to reduce computation time, but scipy's algorithm is capable of solving much larger problems, although required time will grow proportionally."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now define the function which returns a vector of residuals. We use numpy vectorized computations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rotate(points, rot_vecs):\n",
    "    \"\"\"Rotate points by given rotation vectors.\n",
    "    \n",
    "    Rodrigues' rotation formula is used.\n",
    "    \"\"\"\n",
    "    theta = np.linalg.norm(rot_vecs, axis=1)[:, np.newaxis] #np.newaxis converts this into a column vector.\n",
    "    with np.errstate(invalid='ignore'):\n",
    "        v = rot_vecs / theta\n",
    "        v = np.nan_to_num(v)\n",
    "    dot = np.sum(points * v, axis=1)[:, np.newaxis]\n",
    "    cos_theta = np.cos(theta)\n",
    "    sin_theta = np.sin(theta)\n",
    "\n",
    "    return (cos_theta * points) + ((1 - cos_theta) * v * dot) + (sin_theta * np.cross(v, points))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A short review on Transformations\n",
    "\n",
    "Rodrigues Formula: $$\\mathbf{R}=\\cos \\theta \\mathbf{I}+(1-\\cos \\theta) \\mathbf{n n}^{\\mathrm{T}}+\\sin \\theta \\mathbf{n}^{\\wedge}$$\n",
    "If described by a rotation vector, assuming that the rotation axis is a unit length vector $\\mathbf{n}$ and the angle is $\\theta$, then the vector $\\theta \\mathbf{n}$ can also describe this rotation. Here, rot_vecs = $\\theta \\mathbf{n}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# points_pro = points_3d[:, :2] / points_3d[:, 2, np.newaxis]\n",
    "# points_pro.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def project(points, camera_params):\n",
    "    \"\"\"Convert 3-D points to 2-D by projecting onto images.\"\"\"\n",
    "    points_proj = rotate(points, camera_params[:, :3])\n",
    "    points_proj += camera_params[:, 3:6]\n",
    "    # TODO-for-students: Transformations - project() function:\n",
    "    # In the project() function, would it make any difference if I do translate first, then rotate? Why/why not?\n",
    "    points_proj = -points_proj[:, :2] / points_proj[:, 2, np.newaxis]\n",
    "    f = camera_params[:, 6]\n",
    "    k1 = camera_params[:, 7]\n",
    "    k2 = camera_params[:, 8]\n",
    "    n = np.sum(points_proj**2, axis=1) #n = (r_c)^2\n",
    "    r = 1 + k1 * n + k2 * n**2\n",
    "    points_proj *= (r * f)[:, np.newaxis]\n",
    "    return points_proj"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A short review on camera modelling & radial distortion\n",
    "\n",
    "\n",
    "\n",
    "- Each pixel moves radially away from (barrel) or towards (pincushion) the image center (c).\n",
    "- As a function of distance from $c: r_{c}^{2}=x_{c}^{2}+y_{c}^{2}$.\n",
    "- In code, $n = r_{c}^{2}$\n",
    "- The shift $\\gamma$ can be modelled as: $\\gamma=1+k_{1} r_{c}^{2}+k_{2} r_{c}^{4}$ where ${k}_{1}$ and ${k}_{2}$ are radial distortion parameters.\n",
    "- The modified co-ordinates are:\n",
    "\n",
    "$$\\begin{array}{l}\n",
    "\\hat{x}_{c}=\\gamma x_{c} \\\\\n",
    "\\hat{y}_{c}=\\gamma y_{c}\n",
    "\\end{array} \n",
    "$$\n",
    "\n",
    "- **This is applied before the focal-length multiplier and center shift are applied**: Meaning before $K$ matrix is even applied. But how do we exactly do that?\n",
    "\n",
    "    $$\\mathbf{K}=\\left[\\begin{array}{ccc}\\alpha_{x} & 0 & x_{0} \\\\0 & \\alpha_{y} & y_{0} \\\\0 & 0 & 1\\end{array}\\right] ; \\qquad      \\lambda {p} = \\mathrm{x} =K[R \\quad t] \\mathrm{X}$$\n",
    "\n",
    "    $$x_{final} = \\gamma \\left(\\frac{f_0X}{Z}+c_x \\right)\n",
    "     \\qquad \\color{red} \\bigotimes \\textbf{wrong}$$\n",
    "\n",
    "    $$x_{final} =  \\left(f_0 \\left(\\gamma\\frac{X}{Z} \\right)+c_x \\right)\n",
    "     \\qquad \\color{surd} \\checkmark \\textbf{correct}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summing it up\n",
    "Let $\\pmb{P} = (X, Y, Z)^T$ - a radius-vector of a point, $\\pmb{R}$ - a rotation matrix of a camera, $\\pmb{t}$ - a translation vector of a camera, $f$ - its focal distance, $k_1, k_2$ - its distortion parameters. Then the reprojecting is done as follows:\n",
    "\n",
    "\\begin{align}\n",
    "\\pmb{Q} = \\pmb{R} \\pmb{P} + \\pmb{t} \\\\\n",
    "\\pmb{q} = -\\begin{pmatrix} Q_x / Q_z \\\\ Q_y / Q_z \\end{pmatrix} \\\\\n",
    "\\pmb{p} = f (1 + k_1 \\lVert \\pmb{q} \\rVert^2 + k_2 \\lVert \\pmb{q} \\rVert^4) \\pmb{q}\n",
    "\\end{align}\n",
    "The resulting vector $\\pmb{p}=(x, y)^T$ contains image coordinates of the original point.\n",
    "![radial_distortion_1.png](./misc/radial_distortion_1.png) \n",
    "![radial_distortion_2.png](./misc/radial_distortion_2.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fun(params, n_cameras, n_points, camera_indices, point_indices, points_2d):\n",
    "    \"\"\"Compute residuals.\n",
    "    \n",
    "    `params` contains camera parameters and 3-D coordinates.\n",
    "    \"\"\"\n",
    "    camera_params = params[:n_cameras * 9].reshape((n_cameras, 9))\n",
    "    points_3d = params[n_cameras * 9:].reshape((n_points, 3))\n",
    "    points_proj = project(points_3d[point_indices], camera_params[camera_indices])\n",
    "    return (points_proj - points_2d).ravel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A short review on Structure from Motion\n",
    "### Residual\n",
    "In our lecture, in the residual vector, we nicely wrote the elements in order: 11, 12, 13.., 1N, then 21, 22.. and so on till MN. However, notice that it is not the case here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# print(points_2d[:5], camera_indices[:5], point_indices[:5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "M -> camera, N -> 3D point (in our lectures, NOT in this code)\n",
    "![sfm_residual_1.png](./misc/sfm_residual_1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that computing Jacobian of `fun` is cumbersome, thus we will rely on the finite difference approximation. To make this process time feasible we provide Jacobian sparsity structure (i. e. mark elements which are known to be non-zero):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![sfm_jac_2.png](./misc/sfm_jac_2.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.sparse import lil_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bundle_adjustment_sparsity(n_cameras, n_points, camera_indices, point_indices):\n",
    "    m = camera_indices.size * 2 \n",
    "    #TODO-for-students: Jacobian - bundle_adjustment_sparsity() function:\n",
    "    # m above is not \"M*N\" (*2) unlike our lecture notes. Why is that so?\n",
    "    \n",
    "    n = n_cameras * 9 + n_points * 3 \n",
    "    #TODO-for-students: Jacobian & Parameters - bundle_adjustment_sparsity() function:\n",
    "    # 1. Why are we doing `n_cameras * 9` here instead of `n_cameras * 12`?  Recollect :Every\n",
    "    # individual motion Jacobian was (1*)12 in our lecture notes.\n",
    "    # 2. Ignoring the scale parameters, what was the number of unknown parameters\n",
    "    # in our lecture notes in terms of `n_cameras` and `n_points`? What is it here in \n",
    "    # the code? Is it different? If so, what is and why? Link of notes : https://www.notion.so/saishubodh/Stereo-Structure-from-Motion-9fdd81e4194f4803ac9ba7552df56470#6208437ceaaa47f49cfa0ad736de982a\n",
    "    \n",
    "    \n",
    "    A = lil_matrix((m, n), dtype=int)\n",
    "\n",
    "    i = np.arange(camera_indices.size)\n",
    "    for s in range(9):\n",
    "        A[2 * i, camera_indices * 9 + s] = 1\n",
    "        A[2 * i + 1, camera_indices * 9 + s] = 1\n",
    "\n",
    "    for s in range(3):\n",
    "        A[2 * i, n_cameras * 9 + point_indices * 3 + s] = 1\n",
    "        A[2 * i + 1, n_cameras * 9 + point_indices * 3 + s] = 1\n",
    "    \n",
    "    # TODO-for-students: Sparsity, Residual Vector & Jacobian - bundle_adjustment_sparsity() function:\n",
    "    # Explain what you understand from above 6 lines of code by coding a \n",
    "    # simple toy example yourself to illustrate how it is different from what you've learnt in class.\n",
    "    # ([Coding toy example + elaborating in words]- both are compulsory.)\n",
    "        \n",
    "    return A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### THAT'S IT! Now we are ready to use inbuilt library functions!\n",
    "Now we are ready to run optimization. Let's visualize residuals evaluated with the initial parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = np.hstack((camera_params.ravel(), points_3d.ravel()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "f0 = fun(x0, n_cameras, n_points, camera_indices, point_indices, points_2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f2a61098fd0>]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAx6UlEQVR4nO2dd3hUVfrHv296QgkkEFqA0JEmQqQLSEdAXCsqyq6uqGtddV2wrOuqK5Zdy+KiLO5v7ajYEASkCigtAUMPJQQINQQSQglp5/fH3Elmkjsz987cPu/nefJk5t4797xz597vOec973kPCSHAMAzDOJMIsw1gGIZh9INFnmEYxsGwyDMMwzgYFnmGYRgHwyLPMAzjYKLMNsCTRo0aibS0NLPNYBiGsRWZmZmnhBCN5fZZSuTT0tKQkZFhthkMwzC2gogO+trH7hqGYRgHwyLPMAzjYFjkGYZhHAyLPMMwjINhkWcYhnEwLPIMwzAOhkWeYRjGwbDIM45l8fZjKDh3yWwzGMZUWOQZR1J4oRT3fbwZd33Ak+uY8IZFnnEkZRWuxXCOnLlgsiUMYy4s8gzDMA6GRV6G/OJLGPb6KhwsOG+2KUyQHC8qAQCcOldqsiUMYy4s8jLMzzqKnFPn8b9fcs02hQmSDQcKzDaBYSwBizwT1hSXlOHFBTtxqbzCbFMYRhdY5Jmw5q1lezFn7QF8semw2aYwjC6wyDO2o6SsAne8vwHZx4sBAGv3nkLatIXYfqRI9bnKKioBABWVQlMbGcYqaCbyRBRJRFuIaIH0PomIlhLRXul/Q63KYsKbzYfOYM3eU3hu/nYAwLJdJwAAGw+cNtMshrEkWrbkHwGwy+P9NADLhRAdACyX3tsCIbhVF27wL844FU1EnohSAYwDMMdj80QAH0ivPwBwnRZlBcv1//4ZX2/OU/UZAik6rryiEk98mYWc/HPBmMboAJGy307pcQxjV7Rqyb8J4EkAlR7bmgghjgGA9D9F7oNENJWIMogoIz8/XyNzarP5UCEe+yJLl3NvO1KEeZl5+KNO52fkcXe48qX8NNwaZ5jahCzyRDQewEkhRGYwnxdCzBZCpAsh0hs3ll1s3PKYKS57ThTj/bUHTLTAeGr2sBZuPQYA2Hn0rN/PVVYKnDxboptd/jhZXILf/t9GFF0sQ/bxYtPsYMIPLVryAwFcS0S5AOYCGEZEHwM4QUTNAED6f1KDsiyNGR3/8W+vxQsLdppQsnlsO1IIADhY4J2XptJjLEXut3hj2R70+fvyqtmwnug9DDNr1X6sys7HvMw8jH5zNfrPWKFvgQwjEbLICyGmCyFShRBpACYBWCGEmAxgPoAp0mFTAHwXallWxcxx2tKKysAHOYyvMo8AAI6rbA2vzHa1M/KLzU8/zCGbjFHoGSc/A8BIItoLYKT03tHwGJ7+VFQKZJ8oDukcgr33TBgRpeXJhBCrAKySXhcAGK7l+a0Li4ZRlJT5Tj8ghMBH6w9iVJcmshWu0mgphnESPOPVD9wytxdHC0vw7LfbcfcHm7y2V1QKrNh9AsckX/ym3DO1PsvVNONUWOQ1hOsE/fFX8ZZXusYnCi+UeW2fvToHd/0vA6ekUMtPNxxUdD6GcQIs8hpghQmyp8+X4mKptTMppk1biLRpC6ver88pQLkOA8dCeFe4h3l1KCaMYZHXEDNnT/Z6YSmue+dn08pXwy/7XAnFJs1ejzeX7VX12Q05oeenkauTX1iws2pGdN6ZC/hmi7rZ0QxjVVjkHUSoUSdGMWPx7qrXM1fuU/XZeX5SU7jF+0jhRew+7v9a/OGTTPxndY73uTNd575h1i/44+dZqOQwR8YBsMjLYAX3CxMacwPkh/9h23G89MMu2YlRJy0QR68nH60/iBtm/WK2GYxBaBpC6TSUOl+4TvBPZaXAfo2St2niEPP4wRZtP67FGW3Fs99uN9sExkC4Je/B/KyjSJu2EHnSQN2RwouqPs+BGtVk5J7G8H+swsXSCsxcuQ8j31ite5lKe2DlPtwwv+z3XheWK2/GCbDIe/DtFtd0ebdve9H241i791TAz7F7pzYvLNyF/fnnsfv4WWw5VDsuPVj8DW772lXz9zl02ne0zcXSCtMr69LySuw+7j/ZGsMohUU+AGoeNo65VkYolWI4XOK/fr8DY95cg2NF6nqSgTh17hK6PbdE03My1ifsRL7oYhn+u/aApqs/8UpSvnl47haszNZunQAjK1KzftdMaUZu0cWyAEeqI/3FZTh3qVzTczLWJ+xE/qlvtuFvC3bKT20P8aEO99wo248U1XJvHT6tbWvUH75/PnW/q1VWi+K2A6MFYRddU3ihFIDL76kEJQ88P4suxv9rLQAgd8Y4v8eFoqF6y68V9N0KNjDOIexa8rrCD2c1fpqhWrshqorU5awMY2/CTuT9dYFZJIyh5opOVsXs+4HdNYwWhJ3Iu/HXJfb0ratpnPM0eA9M8DloVSJ3yBgnEbYiryXrc1yTaDIOahcPzmiHXVvEvIIVowVhJ/JKHni1D9eJs87OdWIljIx8MatyMOo7Kg0+YOxN2Im8G7nHSO6hVva8cYsrFNSErvo7VotfYe3eU2ET3ZIdIFMn4wzCQuQ9hUFJK11tvLtd3QG6YqGLomYC0OzVOSircNleqsOCJkc98iGt2H3C77EWuoSMjQkLkZdFoY4rOYwfxtDQ+/ot2HpM8bGejYC5Gw9pbsv7aw9Uvb7rfxk4W1I7nNSojsTWI4VIm7YQO44WGVQiYwbhK/IysFZriE4+j+IS46bl69GSr0lFhXl33ZIdrp7EuLfX4osM+fz78zLzsGDrUSPNYjQm7ERej1YjR0GEhpqr5/dYh3SpjBoT2HSgeinFJ+dtlT3miS+z8OCnW4wxiNGFsBB52QFVP51iT9FWlNbAGdqiKYcKzpttQsgY8buaeetcLLP2wu9O4u7/bcKrHsteGklYiDxjPGcu6JO6wF+Vq4VgWqnCtpItwSCEQPe/LsGnG7Qf27Aby3efxL9X7Tel7LATeX/PjVx43tfSQiJG8NayvfgiwNqkTkRVCKWOdtQqyySVdUoIZ3mlQHFJOf7yHS83aCZhIfJyj6rSBynrcGHAY4Z2SgEApDaMV26TEHhyXhYyD1b7Rd9YtgdPfiXvG2Vc+BNeLbTRaIH1H/dv36b8saKLqOA0H5YgLEQ+VLKPFyNt2kJs9rGMXZ3YSABA28Z1FZ/zYlkFvsjIw+Q5GxUdv3zXCeSeCuzn/u7XIxg4Y4WtHjArWeqpuabNeLV59pzCC6Xo//IKXPP2GrNNsRzHi0oMLzPsRL5cCotTM+g0+k3XItTfybhuzpaUGSIGd3+QgaGvrwp43LSvtuFI4UVcKnd9vw05BTgfJqsBaf0zeJ7vu1+P4OUfdmlcgjNxh7nm5Nt/8F1rftpz0vAyw27RkM2HCgEAby/fi6slN0tNlIr22ZIy9Pjrj7isWf2g7dGrS/7jjhNYuO0Ylu48gRGXpaBxvVhsOHAaKx4fCsC13mdCTCQSYsy/BTyv91vL9qK8shKPj+oU0nmCtsXH9kfm/goAmH7NZaEXUqO8gwXn0bJhAiIivFvwdh94ZaxBWLTk5fyeJWXqJ7rUPEuRFEGy65jyxb7d6N0lf/TzX7F0p2uyy65jxfhs42Hk5J/H377fCcC13uf4t9fqaoNSFm2vnpH6xrI9+NeKfSZaU42cyGrtBttzvBhDXluFd1dXR15sO+KagWq39VjLKiqxPqcAV7260mxTLMHDn23Bkh3HzTYjPEReL2qKQDCyrXVrTa5n4FnJ/ffn6mn1OafO461le/HHz3/V1giVfKMigknvDI1yeY48B9+vn/WLpuXlnXHlssmQWXP4mAn+22ARQqDD04swafZ6s01RhBBC9+ip+VlHce9HmV7bVtdYA9kIQhZ5ImpJRCuJaBcR7SCiR6TtSUS0lIj2Sv8bhm6udijNQqlmv5u1e0/h1Dn/6Yf10iq5Hoo/s99YtkeVyOrBxVJtJuXo5fqa+M7PVa+VRFuFI+7KKlj2nCjGmfOlGlnjn0vlFWgz/Qe0mf6DLucvKavwuqdPe3yvhSryKGmFFg7ZcgCPCyE2E1E9AJlEtBTAbwEsF0LMIKJpAKYB+LMG5VkaIQQmv78B7VPqYtljQwIfb4BNlRZ37paoGAT3G3Koo09eL+wcJumJv2uvpEEz6o3VaFI/FhueGuH3uMOnXUtHni8tR+emysfCLpVXoLS8Ek98mYXICH17g+kvLvNytd39wSZdywtEyCIvhDgG4Jj0upiIdgFoAWAigKHSYR8AWAWTRF6rx6jmA+nvAd138pxGpYZOzQcwbdrCWsecLC5BSr04gyzyRqvf54JGPQI3hqQ1kMqwd9CkNgRafGdbXhEmzPQeR9ry7Eg0rBMT8Nydnlkckm1qqDmW4iv0Oft4Ma6duRYrnhiKFg2Uz7FRi6Y+eSJKA3AFgA0AmkgVgLsikA1lIaKpRJRBRBn5+flamhPAVt/7lDzbQgjYJRT9ZHHglav6vLTctBmegYr9PuuozzkKnhxQMI9AjS1GXA1/ZdhJ+HNO6d+okSvjaFFobiJflFVUVoVb68WnGw7iUnklluo8OKuZyBNRXQBfAXhUCKE43EQIMVsIkS6ESG/cuLFW5gRFVctcwdN9z4eZuFpB3LqCQmU5WVxSFetuFPvztXlQ1aZmqBvrv0P50GdbcP2/XQOeRi7/pyW+KrKqlrw9v1YVRuSnUdsGKS2vrFp/WQ2ZB0+jw9OL0O/lFYo/k3W4EDe/t072mTX7ntVE5IkoGi6B/0QI8bW0+QQRNZP2NwNg/CwACV/L+r22ZLes68JXE8rzPMt21V7V56c98j2RbXlF2JbnvTBDoN+9z0vL8cAnm/0fpDFTP8wMfJACAqVm6JGa6PU+KtI6CmeTzlnY8PqSbLyz0hVSK+ceXesnWuWVxbsxafb6Ws9eIG6YtQ4AAgZPePL0t9uw8cBp7DmuvKFUatBaAlpE1xCA9wHsEkL802PXfABTpNdTAHwXallaQiC8szK4rHB7T6hbG3PCzLW1fIlu/Pn1l+0yrV7Ula01HjorjQuf8nRtGWCYUwZe/c0fKAtBzGau3IfXlmQDkI8ce3mR7/S9e6VxsYLzysVaD077iBr6TFp5bJNM+KyWaBFdMxDAHQC2EdGv0ranAMwA8AUR3Q3gEICbNChLVwI90+XSzfrQZ6EvouA5GerUuUteiyp7RpucOFuCkx4DUtvyitBdagmb5UM3E72/8/nS6kGzvQYMnld/Hev0ZoLhl/3q3SJqKVDRstaK0vJKxET5bwtXVgpsP6J+QqSbMxf0DR0NuSUvhFgrhCAhRA8hRE/p7wchRIEQYrgQooP0/3Tgs2nLwq3H0PuFpSiv1GYA5XNpibTdGq5yLwRwy3vrcPucDVXb/vFjdtXrIa+t9OoFTJi5FrfPWY9L5RWat/TNqjKs2ppdtN382YpMNXKdBSXRkErvrk25p2u5bzs+syjg55QEBZiJ+YlLdOS5+TtQcL5Utrvk1yeuoeb4igHfLi2eXF4psL9GIqf/rKmelSrXRf15XwGG/+OnkCeg2BG9qwOjO0fWrN6sidx8D38RblWPuMKLfNO761TbVNOGE2fVz1LW+55zeFqD4K6eli3LMh9hWJsPhlb76yHwvmy1EquyjQuzNYJnv3UtqGH36BqtyMg9jdmr5cfK1M6DcF9TvXuKWzxa8r//MEPXsoLBkS35gnOXcOt/1lfdFEpDmJTUqCdV1NSfbjiEp77ZJrtvZxBJzfRGbcWxePsxfJ91DO/c3iukcq00tFB0UZ9lCxnfePZ2b5Ra0wkxUeja3HtG6+zVOarOq6bePFRwIeAxWYcL0bFJPcTHRHpt9zf4qwS9KyFHivzIN1b7HNF243kDLNt5Ak0T4zz2+b49+vx9uWI75AT+SOFFnL1YhtJy67eaAVf3c3/+OQxo1wgAsPv4WRwquIBRXZvivo9dIZ5Nvt+Jv0zoonnZS3eewKHTgR8+PZm5Yq+p5d/zYQbW7S/A9udHm2qHnnR+tvZs1Ge+NXbJwMGv+c+ceercJUx852dc070p/n17b03L1ruR4ziRP+PDB+8Pdxerf9tkAKHVrMP/8ZPf/QNnuCZYjO/RLOgyjGTsW2tw+nwpcmeMw3PfbccH6w7WOua/Px/QReTvsUDX9/Uf9xhSjq9mhTtdtNWxsrtJCxF1JxyrGf5rBxwj8kUXyvDqkt1YIzM5Qvb+k7krteg25Wgwtd5KuCtMIYSswGuB3EOoNFRy1qrg5jpYDSuLpBKsaL7bTRuqyH+9OQ9XpiVpYJE8ensrHSPyl//tR1XHy92UGkVaKkLrZFp6sy6I6eGh4F6JKRCvLA7NH2oV7L6uqxVxX9FQRXT26pwqkbfS+JFSHB5d40Lpeq4bc/UN5d95tHqw1ei8NKESyqzFYJifddTQ8sxmy+EzeHnRLqzcLT/3wYyJQEag5+pXy31cS8vBIZShE8hPXhO9amvP1eutGsXhKxeInu1Mo3sJRlNRKfC/X3L9HnPi7CW891MOfve/6tzjj3+ZVfX6Bo1XpLICJWUV6PbcEl3OvWZvdahtqGM7Wi/5WBO9o2vCQuTl8OcDNWJtzVCmQevJ5Pc3yG63YS/VNE6fL8WYN1fjYIFrfEaLXluughA/u3EpiHWWlfKThvMptE5vYXQ6krAVeX9ombbAjmw/UoSCc5e8Kju9PcbTvtoqnxHURhwquICii2VYuO0Ydh8vxuzVOSgpq8DGA4Zn9Ah77DSQzSGUCrDDTE07Mf5frlw5Y7s1rdoW6KGpDLFLO1dlDnorMvi1lUhtGI97h7QD4JruLhcD7kTU/PpnzpdiXU4BBrRL1s0erXO4u093pNB+qUQc0ZL/ZrP6hahtVNGbhpoEXT9sN36BYiviOWt4x1FtY6oLdc5WGApqosXu/yQTf/hkc8Dl/kJB6fMdDplcHSHypUG05M9csObAp12xW0ionrgFRusF1Gf95Iw5Ae6KcPSbq022xJjxt0BkhJjHKhCOEPlgemZarAfKeGBAg+iXfb5XAbIiYdBItCxKc0OZvTSfmyIdG53OEHl2vuhOoGtsRE74ez/SZnlCo2CRNw+tQ5R9VQZZhwtVn+toUe0kh5cq9OsJO0LklSwcwISGFRo8l2wywG6Fa2VljLg+Rv0EavNkAUCunBdBxwaBI0SeHyrzMaLVauXMnaFGFymBe6zW4+P1weVzqjl3Qs+7xxkizze/7gQS2HAfeH3bIyWx+35kb43zCSZ1AqH2im96NpKcIfKs8bqzKECI5OIwXw91vUxqBq3D8/KLnZm/RglKFvXwwuKiUNM8Pce0HDEZyioj5E7GAG+E7XCnOb45PdVr++nzLjHWeub0sl32yC2vB4EW9fBk7sZDQQ2ImomeLXlniLzZBoQBAX3OYfgjuNMcv7J4N65Ma1i1/dQ5fSYtFV0sgxCCGzUBmPa1/JKbVqG8UqC4xDs+X+s5FZ6wu4ZRBDfk/eP5jOp5P6qZhWxFzpaU4fBp+6UG8PxJr525FodPX8Dnmw4Fda47/7uxaoU4N9ySD0AEq7zuhMP071DwvDp63o9ngkxtEKwgackXGYfx5LytZpsRMlvzivDe6v34eL3511QJ3JJnFBHQJx/mdUCmx9T0SB0nbgRb1/75q9BcGOUazFGws8Cv26/vmgccXcOYTkWgu5Ar2ir0bHQE06Panx96PvSvt6hPAugkPBdwAewVtu0IkWd3jQFwS14xWi5YUZNgLrPSldH8rYD0+abDOFZ0EZkHT2NbnvrsmrNXOyO5mhutJUfPEEpHiDxrvP4Eugmz8gqNMcQi7POzWpCei86obciria1v99QPXu+/ysyrep158AwmzV6PG2atw4SZa1X3KF5ZnK3q+GC4VF6BLYf0zejoRmvJYXdNAOzUdbIrhQGy5IXbuOyIf6pbN1gr1Irru0GkJ37vp/3Yd/IcXv/RW5iPeyTW+myja5EXpZWInuukZh48jcOnL+DZb7fjN/82Zi1crcNY9Xx8HBFdY0QGxHDnlwADT+WV1s0rozUzPVIYGE3NO33Wqv34MuMw5t0/AEl1YlBw7hK+2XIEdw9qgzeX7cXKbHXT7i+VV+DlRbvx8qLdfst+6pttaJ2cgNvnbMCrN/RQ/0U05IZZ64L+rBACZ0vUZ6wMtDC7WlbuPok2g9poek43jhB5xnzCaUbs6z/uMa3sb389it9c0QLbj5zFupxTeGelq6Xe64WlmP/gQLy2JBtr9p5C8wbxeGu5+srozHnfglczf9Htc1yLvn+/9ajqcqzCf1bn4O0V+8w2A5mHzuAu6CPyZKX45/T0dJGRkaH6c3M3HrL8LDeGsQPtU+r6HW9g9CP7xTGIjYoM6rNElCmESJfbp7tPnojGEFE2Ee0joml6lHHqXPgmbmIYLWGBN49Oz+iz6LuuIk9EkQDeATAWQBcAtxJRF63LaZmUoPUpGYZhHIHeLfk+APYJIXKEEKUA5gKYqHUh9eJ4aIFhGEYOvUW+BYDDHu/zpG1VENFUIsogooz8/OAmkSTGxwRvIcMwjIPRW+Tlgkm9RnqFELOFEOlCiPTGjRsHVUgUL/LKMAwji94inwegpcf7VAD2jbdiGIaxGXqL/CYAHYioDRHFAJgEYL7OZdaiFQ/MMgxjcf5zp2wEZMjoKvJCiHIADwJYAmAXgC+EEDv0LBMAbuqdipy/X4MGCdEAgCb1Y/UukmFMI7mOa0xq7tR+JlvCBMtvB6RhZJcmupxb9zh5IcQPQoiOQoh2QoiX9C4PAPq3S0ZEBCE2KkKywYhSGUZ/Nj41vOr1p7/vi4eGtcfP04Yhd8Y49GubjAUPDULDhGhs/euooM5/dSff42JN68cFdU4mMCk6NkQdkaDMFyzujNNIqR+Heff1x8zbrsCA9o3w+KhOiIuuniXZrUUitvxlFOrHRWP3C2MwvkczVeefM+VKn/s4R5R+3N6ntW7ndqTIc+phxsmkpyVhfI/mAY+Li45UHXkWGUG4s79+gsPIQzoqsSNFfmD7Rqo/07lpPR0sYRhzCSYl7oTL5SsQz57xi9d1Q4/UxGDNYgzEcSKf/eIYpNRz+Q6Vdi7/NLoT7h/aTj+jGMYkgunUprdu6Hd/w4RoTO7XGnOn9sOL13ULzjDGMBwn8sFkcXvg6vZVg7QM4yiCUHl3679r8/pe2wWAXX8bg3XTXYO/CTFRmNyPXTtWh5VNggdpGaaaNU9ejXn3Dai1PT4m0mugNxyx2wz7sBB51m+GUUfLpATEx3iLOTeEXOz7+zXY9PQITc+pZ7XhaJFXc1MGOrR1cgJGXKbPZAWGYeyFnVrzjhL52Xf0DvqzgWbFPjSsA9qn1K16f62PCASG0YtJV7YMfFANtF3kPvSmvFNm5UawyJtDVGTwFz4+Wl1O+jqxnMOeMZZhnVNUf8Zqc0b6tU022wRN4Ja8RXj1xu7o3LQeGko5bEJhQDvvm/PKNP9hZgxjBTRtx7NPvoo6sVF49cYeZpuhCEeL/LDOTbD40cGIVFDrBmrxNG8Q7/V+VNemoZjGMIzNuTldvfvMDBwt8m70aIHUtbm7pmbPhLE+wbgjtXQrckPenjhC5Js3iMfz13ZFhxR1qQk8s+pZzXepN+H2fYPFSmmqYyLVx6df1sw66ToCuU1jLDoh0TPgQgtSG8YHPkhDrHlVVdK4XiymDEhDSx+Lg/hqgXz/0CB8db9rwke8ggkeDTTw7YcL8+7rb7YJmrD2z8PMNsEyiBC7xAHz6Fi0qzC5bytNz/fnMZ2x5dmRmp7TH44Q+WBpXC8WvaU8Ha2T6wQ83kmNX7nQuhYNtGthBNtT6N7CWkmvoiOt84jYvfcVaGhMr1TGl6cm+ky6poTLmtUPfJAKiIz9La1zB+uI0gZIs0RnLorwmyta1Nomd5PdO6RtUOdf8+TVtbYF2+h7YnQn1Iuz93iHP+oF8JE7ORtqMBkxteD2vq3xwsSu6NYiOLHu2zZZU1dShMHXISxEXildm8u3Imfd3stgS7TljVt6Kj729ZsuxxWtGqg6f1x0JLKeG4U6HtPgQ2mT/Wl0pxA+bW0CVWByFbIbvaXhlRu6+90fajs7ULBCqAES/7z5ctntQzs1RoOEGPxxRMegz10/Tr2r9nofvyWhdk9azwowTES+9t3jXhdTCWO7q1tdx87c2DsV3/xhIHJnjMNX9yvzqxMBifHRePVG+YdMLU0cvMxcoIe5YxPzWvJ6xsE/PrIjPryrj//yA5zD3zObO2Mcru+VKrsvRbqfQvl+V3VQv0bFLT5mKBvdoQkTka+NFi6B//5Wn9XV9eTVG1wTOOTEpuZD0Lt1kuw5fA1Aj+1WPXcglAdqVJcmeGtST69tgQTihl6puLWPPeKW/XG1v1mtJvvkQ/lNHxrewWdghJtOHhVcTRdg7oxxpgY+BCPMfdsm+2gokaG/ZdiKvHyLqvZd/OQY366DYZ0DJyzrkyYvlEZR0/7GUkhgKPfY+unDvSo497mCzecxfWxnLxcREWFiT++u7uCOtReYvqXWZBSbj0xKdJBC9t65rRfuHtSmanu9WH1FruY1N5rY6Ai8Nakn7hrYRrZCUFLHtExSFzzwpkJXZmSQzW+5hhK35HUg2BbI5H6t8Ieh7aveK5k5W5M7B1QvqtDJ4K547oxxXvYD8PukyIXIeebo2PDUcMx/cCDioiMDVnD14+V7Sv+69Qqv9w8Na497h7TDGzf3xNhuTdGvrfJK8RULTSt//tqumvfsOjSpizs8FuXorvNyezVTC9ck1BBKJUzs2QJ/mdBF8fEjLmuCDzx6eYsfGezzWDlX2HVXtMD//da1ePmmp0dgkI+lQ6eN7azYpkAQOLpGN/p7JEcK5hoHswqO5wDLnCn+ReDewfLRLUNkWrGefK4is587TE3uJpN7hH99bhTWTx+OHc+PRpP6ceiR2qBq3xgptUOszByDzk3rY3I/7/ji+OhITLi8udcArZu0RnUwa3Jv2ZW9fF0XAFj9p6uRVCcGU/0cEyqerWkAyPrLqFrH3HJlS0U9O7W4fye1LVQ3ft0/GvHMuMt0LwOQf2ZbJsV7PR91YqPw+k2X49Pf9611bKvkBNl1aa/unILcGePQuJ7viW/JdbWbFGd0lFFYiLxbvLweVgXXueYIeFx0JKaN7YxFj1ylqNwXJnb1Pp/H6eY/OLDW8U+O6YykOjF4+prqhyZ3xjivloocfVVk9nNn20yuo+ymrRsbhaaJcbLT49+4pScWP3qVz6iJIR21EZjp1/gWkVbJCdj87Eh00jH0cLRUmY24LAXLHhuMRBnfsBarJckN7rnD7YJtRLvXO9YCXybc2T9NszL8oWQcCXAFDwzw0SJvFWBcIMWP0H8+tZ9f961SCMY6FsNC5N143iPBXuT7hrTzOzli8aPVFcBN6S29ymzRIB5/HNERq/90NXqkNsC3D3gLfWQEYfOzI3GPTKv0+Wu9K4xArh9fkQj92ibhpd90w/M1KiBAvZDEx0Sic1Pf1yIx3iWG7jzovlw4vvj+wUG13Dv+kGsgaZGWoE+bJGQ9NwpzplyJ9ipTZyilVVICPrq7duvT7UIxM+omEGpjyBvVaBW7XSE179mPZa5HqMy4wb+L74XruiHBh9uqb9vk2u7PIDDaJ+/cWSceePoS2zaqg5xT53Frn9pTlYNpLV2Z1hCbcs/g8tRENK4Xi05N6qF5YhwullUgLjrSqwVMRHhkRIeq9z1bNkDbxnWQk3++1nnfndwbFZXVBt3RrzWiIyMw/vJmOFdSjuYN4rFo2zF8s+WI1+fen5KOnPzztVwMnjbc3teYxZf7tEnCe3f0xtBOjdGleX1Zt9NoP9k8u6cm6u6H9odny9pdYbl5eFh7vL1iHwBXeKBeNKobi8/u6WfqdahCI5d8nzYN8cO241Xv7x3cFvXiojCuRqjyoBo9mxiVs48XP3oVso8Xe20LFKtfJzYKN/ZOxYfrDqoqSw2uGa/GKX3YteTdKYP9de+bq5j5OlDqFg7plII5U64EEeHnacOwWcpNMbhDIwxsn+wz9/SKx4fKbh/TrSnG9ai+6SMiCLf1bYX6cdFV32Fs92aYfae3n3/4ZU1wz+C2QUW66DGsNrprU8RGReLO/mlVqSPSGrn+7/zbaHRTkcbA6Nmg/sZCHhvVqWqSnOe9tOqJoUGV5e+Z798u2bSsp+/c1qtWOGuovHx9Dy9Xprvh0SDB/9yVx1RWpp2b1jc9YkgObVfrCkxYtOS7tUjEyux8xf7Jdil1cbSoRFG3anDHxnhz2V4M6Vjd6vCspYkIn/ze/8BoUp0YnD5fqsi2YOnZsgEKL+hbhlI+vKsPsvIKkRCj7vZb8NAgVAqg4zOLdLLMxYTLm+P7rKMBp5+P7d4MKx4fgraNq7MUpjWqU+v3bJ4Yh6NFJbrZqwW+XC7jejTD2ZIyANo1AhLjo3FH/zT0aZOs6p50SroLkgmT13OlKWdctQA8OqIjRnZpoqrVqJRerRoid8a4kM6x8omhuFhaoZFF8tT0/5tJct3YoCJRogJ017V4TPq2SVI1G9pT4OXs+PaBgUhtGI/0F5dpYF1wLH70Kox5c03Qn9dLfpQMlr9wXTe0bRQ4eaCdiCBCXHQk4qIjMG1MZwzplKLJwL0vwkLkIyPIK/TPaiTGR9fy+ZqBEXHQeiIbFmriV1rz5NUBZ3m6UetvVoN73YTHR3bEP5buCfo8ZtwfnvME5Eq34z1L5NKk3S+MNaS8sPLJA/7TmdrvdrEe707uhR8eVhZiqjW+fJ1ys2V94W8gWLEdkhlqWmdJKnoPammQEIN9L43Fg8N8R4b0bNnA5z6zskc6FaN98mEn8m6UXGin39qeA7taMaZbM3Rprm3+baUM9BEbLTf5yhdyUVdqefvWKzCofSNVwq23jkZFRvgVa39ZFt3u4qS6+lVEwWLHhpmtQiiJ6DUAEwCUAtgP4HdCiEJp33QAdwOoAPCwEGJJaKYyWvPObb2wcOtCs83QDC0Gr7R4AAe0a4QB7dRlLayjchDaSBJiovDKDd0xqIPyHhHgSoY3/LIUREVYqy05d2o/TJq93ud+vTXY6MZjqFd/KYBuQogeAPYAmA4ARNQFwCQAXQGMAfBvItJvZIFhEJxAp0srg5mNUeub5s4YV7UamhpuubKV35XDfi8zL+PmK1siuW6s7AxhM+kXYIb4HR4zeIOdTOf3c3ZKUCaE+FEIUS69XQ/AndB5IoC5QohLQogDAPYB8D83nwmZdo2Dj0Lo3LSeJWOKtcA9Q3nmbVfg/qHtTLamGs90FZ4VlHsCnV4rCH18d1+snz4cP0+rXr92fIiuu2fGd6lK9KUXSgexlfDWpJ749B75GbWBUh+Eip3j5O8C8Ln0ugVcou8mT9pWCyKaCmAqALRqpe2CuWoZ1L4RVuw+idbJCViz114DThufHh5Sl3/xo76z99mdB65uj6GdGqNHagOM79Ecs1btN9yGAe2S8cv+AgzrnIKP1vufTfnu5N74ZsuRkCptf8THRFalS4iLjkBJWaUmA84pGqSQ8EeLBvHY+ldXcriy8krc93Em7h0SXKVtZoPGV9oEvQioCkS0DIDcHfC0EOI76ZinAZQD+MT9MZnjZcdIhBCzAcwGgPT0dFPHUX43MA0TLm+OBVuPmmlGUASbiOreIW1xrNDaE3WUUk8aPJzYszn+OKIjhr6+CkDtENqf/jQUQ15z7TOqHu/avD5+2V/glU2yp4+w3qaJcZbqcSila/NEzH9wIPKLL2na6vbEc4D4y/sG6FKG5z3RJchFvP211rsaHJgQUOSFECP87SeiKQDGAxguqoNW8wB4ruiQCsDyyklEftONOpHpY41JE2sEfdok4a1JPTGqS9OqWZpytE6ug85N62H38WJc27MFNuWe0d02uXDuxIRozLztCjz46RbdyzcKK89HCYZ/3abt+s4PD2sfcFKf1oRUGhGNAfBnANcKIS547JoPYBIRxRJRGwAdAGwMpSytcGfza1jHWoNBjDrap9SeaQq4uuHxChYUd6+/OaFHs4BLC2qJrxae0X5axjdREYT7h7bDgocGaZ4zKNpggQdC98nPBBALYKnkv14vhLhPCLGDiL4AsBMuN84DQgh95+0rZPrYyzCma1N0be47xYENJ9GFFfPu6y+bTsCTQJL5u4Ft8LuBvjJ1uv7zbRCeEBH+PEa7laDMJiSRF0L4nEInhHgJwEuhnF8PYqIiVC2ywViPdI3XzfX0wcZERuDR4R1x/lI5btNgYpQbXxVGmpSZM5iwRi24tU8r/N/PuYiO5J6EllgpZsO6MzAYxiDcg2tz7kxHRAQhMSEar954uS5l1Xz4u7VIxE9/Gqp72J4vnh3XBdPGdjbcT8wYB4s8E/Yk140NOZNoIPy5AN159s0gIoIQG8HzFJ0MV98ysC+WYZhQ8OXqM0NbuCXvByv51Rh74y/7KeMsDrx8jdkmeMEizzAG4p5FPaxzismWMHphtZnyLPIMYzB6+/8ZxhMWeRniol1DFUbnmGCcyyPDO+BcSTlu7dMy8MGMbUhtGI+8MxfNNsMvLPIy3JzeEgXnSnHPVW3NNoUJkggpt3xDi6S5bZAQg9du0icsk2H8wSIvQ3RkBB4e3sFsM5gQaFQ3Fi/9phv7vpmwh0MoGcdye9/WaJboe6ELhjGCZY8NMbV8FnmGYRgd8ZVMzyhY5BmGYYLEDskMWeQZhmEMwoxKgUWeYRjGwbDIMwzDBEmzxOCW3TQSFnmGYZggSaoTY7YJAWGRZxiG0Yj5Dw4024RasMgzDMMESf147xnVVlzInEWeYRgmSJ6b0MVsEwLCIs8wDBMk9eKskRvJHyzyDMMwDoZFnmEYxsGwyDMMwxiEe60KI2GRZxiGMYjfDkwzvEwWeYZhGIOIjTJ+tTkWeYZhGAfDIs8wDONgWOQZhmEcDIs8wzCMg+GFvBmGYXQmqU4M+rVNMqVsFnmGYRid2fzsSNPKZncNwzCMg9FE5InoCSISRNTIY9t0ItpHRNlENFqLchiGYRh1hOyuIaKWAEYCOOSxrQuASQC6AmgOYBkRdRRCVIRaHsMwDKMcLVrybwB4EoDnOuQTAcwVQlwSQhwAsA9AHw3KYhiGYVQQksgT0bUAjgghsmrsagHgsMf7PGmb3DmmElEGEWXk5+eHYg7DMAxTg4DuGiJaBqCpzK6nATwFYJTcx2S2CZltEELMBjAbANLT02WPYRiGYYIjoMgLIUbIbSei7gDaAMgiIgBIBbCZiPrA1XJv6XF4KoCjIVvLMAxjUe7s39psE2QJ2l0jhNgmhEgRQqQJIdLgEvZeQojjAOYDmEREsUTUBkAHABs1sZhhGMaC/G1iN7NNkEWXyVBCiB1E9AWAnQDKATzAkTUMwzDGo5nIS615z/cvAXhJq/MzDMMw6uEZrwzDMA6GRZ5hGMbBsMgzDMM4GBZ5hmEYB8MizzAM42BY5BmGYRwMizzDMIyDYZFnGIZxMCzyDMMwDoZFnmEYxsGwyDMMwzgYFnmGYRgHwyLPMAzjYFjkGYZhHAyLPMMwjINhkWcYhnEwLPIMwzAORpfl/xiGYcKFWbf3Qlx0pNlm+IRFnmEYJgTGdm9mtgl+YXcNwzCMg2GRZxiGcTAs8gzDMA6GRZ5hGMbBsMgzDMM4GBZ5hmEYB8MizzAM42BY5BmGYRwMCSHMtqEKIsoHcDCEUzQCcEojc4zGzrYDbL/ZsP3mYQXbWwshGsvtsJTIhwoRZQgh0s22IxjsbDvA9psN228eVred3TUMwzAOhkWeYRjGwThN5GebbUAI2Nl2gO03G7bfPCxtu6N88gzDMIw3TmvJMwzDMB6wyDMMwzgYR4g8EY0homwi2kdE00y25b9EdJKItntsSyKipUS0V/rf0GPfdMnubCIa7bG9NxFtk/a9TUQkbY8los+l7RuIKE1D21sS0Uoi2kVEO4joEZvZH0dEG4koS7L/eTvZL50/koi2ENECG9qeK5X7KxFl2ND+BkQ0j4h2S89AfzvZ7xMhhK3/AEQC2A+gLYAYAFkAuphoz2AAvQBs99j2KoBp0utpAF6RXneR7I0F0Eb6HpHSvo0A+gMgAIsAjJW2/wHAu9LrSQA+19D2ZgB6Sa/rAdgj2WgX+wlAXel1NIANAPrZxX7pnI8B+BTAAjvdO9I5cwE0qrHNTvZ/AOD30usYAA3sZL/P72VEIbp+AdfFXOLxfjqA6SbblAZvkc8G0Ex63QxAtpytAJZI36cZgN0e228F8J7nMdLrKLhm2pFO3+M7ACPtaD+ABACbAfS1i/0AUgEsBzAM1SJvC9ulc+aitsjbwn4A9QEcqHk+u9jv788J7poWAA57vM+TtlmJJkKIYwAg/U+RtvuyvYX0uuZ2r88IIcoBFAFI1tpgqSt5BVytYdvYL7k7fgVwEsBSIYSd7H8TwJMAKj222cV2ABAAfiSiTCKaajP72wLIB/B/krtsDhHVsZH9PnGCyJPMNrvEhfqy3d930v37ElFdAF8BeFQIcdbfoT5sMc1+IUSFEKInXK3iPkTUzc/hlrGfiMYDOCmEyFT6ER92mHnvDBRC9AIwFsADRDTYz7FWsz8KLjfrLCHEFQDOw+We8YXV7PeJE0Q+D0BLj/epAI6aZIsvThBRMwCQ/p+UtvuyPU96XXO712eIKApAIoDTWhlKRNFwCfwnQoiv7Wa/GyFEIYBVAMbYxP6BAK4lolwAcwEMI6KPbWI7AEAIcVT6fxLANwD62Mj+PAB5Us8PAObBJfp2sd8nThD5TQA6EFEbIoqBa0Bjvsk21WQ+gCnS6ylw+brd2ydJo+5tAHQAsFHqFhYTUT9pZP7OGp9xn+tGACuE5OQLFams9wHsEkL804b2NyaiBtLreAAjAOy2g/1CiOlCiFQhRBpc9/AKIcRkO9gOAERUh4jquV8DGAVgu13sF0IcB3CYiDpJm4YD2GkX+/2it9PfiD8A18AVCbIfwNMm2/IZgGMAyuCque+Gy++2HMBe6X+Sx/FPS3ZnQxqFl7anw/WQ7AcwE9Wzk+MAfAlgH1yj+G01tH0QXN3HrQB+lf6usZH9PQBskezfDuAv0nZb2O9R9lBUD7zawna4fNpZ0t8O93NoF/ul8/cEkCHdP98CaGgn+339cVoDhmEYB+MEdw3DMAzjAxZ5hmEYB8MizzAM42BY5BmGYRwMizzDMIyDYZFnGIZxMCzyDMMwDub/AafIDlOJfePCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(f0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = bundle_adjustment_sparsity(n_cameras, n_points, camera_indices, point_indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "from scipy.optimize import least_squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Iteration     Total nfev        Cost      Cost reduction    Step norm     Optimality   \n",
      "       0              1         8.5091e+05                                    8.57e+06    \n",
      "       1              3         5.0985e+04      8.00e+05       1.46e+02       1.15e+06    \n",
      "       2              4         1.6077e+04      3.49e+04       2.59e+01       2.43e+05    \n",
      "       3              5         1.4163e+04      1.91e+03       2.86e+02       1.21e+05    \n",
      "       4              7         1.3695e+04      4.67e+02       1.32e+02       2.51e+04    \n",
      "       5              8         1.3481e+04      2.14e+02       2.24e+02       1.54e+04    \n",
      "       6              9         1.3436e+04      4.55e+01       3.18e+02       2.73e+04    \n",
      "       7             10         1.3422e+04      1.37e+01       6.81e+01       2.19e+03    \n",
      "       8             11         1.3418e+04      3.70e+00       1.29e+02       7.96e+03    \n",
      "       9             12         1.3414e+04      4.23e+00       2.64e+01       6.19e+02    \n",
      "      10             13         1.3412e+04      1.88e+00       7.59e+01       2.63e+03    \n",
      "      11             14         1.3410e+04      2.10e+00       1.77e+01       5.01e+02    \n",
      "      12             15         1.3409e+04      1.04e+00       4.04e+01       1.32e+03    \n",
      "`ftol` termination condition is satisfied.\n",
      "Function evaluations 15, initial cost 8.5091e+05, final cost 1.3409e+04, first-order optimality 1.32e+03.\n"
     ]
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "# So far: method='lm'\n",
    "res = least_squares(fun, x0, jac_sparsity=A, verbose=2, x_scale='jac', ftol=1e-4, method='trf',\n",
    "                    args=(n_cameras, n_points, camera_indices, point_indices, points_2d))\n",
    "t1 = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Setting `scaling='jac'` was done to automatically scale the variables and equalize their influence on the cost function (clearly the camera parameters and coordinates of the points are very different entities). This option turned out to be crucial for successfull bundle adjustment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization took 22 seconds\n"
     ]
    }
   ],
   "source": [
    "print(\"Optimization took {0:.0f} seconds\".format(t1 - t0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's plot residuals at the found solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f2a5ba3fdd0>]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvjklEQVR4nO3dd3gU5doG8PtJI6G30MHQEVARIl0QpAqKnqMe9NgL2BU9elBsx3b4/A52j4qKvaAo6odYKCqiCAYJvWOEUEIACYQSSPJ+f2zJ7O7M7szuzpbh/l0XF5vZKe+2Z9553jKilAIRETlTSrwLQERE9mGQJyJyMAZ5IiIHY5AnInIwBnkiIgdLi3cBtBo2bKhycnLiXQwioqSydOnSPUqpbL3nEirI5+TkIC8vL97FICJKKiLyh9FzTNcQETkYgzwRkYMxyBMRORiDPBGRgzHIExE5GIM8EZGDMcgTETkYgzwRUYQOHD2OL5bviHcxdCXUYCgiomR010fLMWdNEU5uUgvtG9eKd3F8sCZPRBShnSVHAABHj1fGuSSBGOSJiByMQZ6IyMEY5ImIIpTIt8pmkCciihKReJcgEIM8EZGDRSXIi8g0EdktIqs0yx4Wke0iku/+d040jkVEROZFqyb/JoAROsufVkp1c/+bHaVjERElFMfn5JVSCwDsi8a+iIgoeuzOyd8iIivc6Zx6eiuIyDgRyRORvOLiYpuLQ0R0YrEzyL8EoC2AbgB2Apiit5JSaqpSKlcplZudrXsfWiIiCpNtQV4pVaSUqlBKVQJ4FUBPu45FRBRPCZySty/Ii0hTzZ8XAFhltC4RkRMkYj/5qMxCKSIfADgLQEMRKQTwEICzRKQbXCe5AgDjo3EsIiIyLypBXil1ic7i16OxbyIiCh9HvBIRRUglcEd5BnkioigRJF5SnkGeiMjBGOSJiByMQZ6IyMEY5ImIoiQR+8kzyBMRORiDPBGRgzHIExE5GIM8EVGEEngsFIM8EVG0sOGViIhiikGeiMjBGOSJiCKkEvi2IQzyRERRwgnKiIgophjkiYgcjEGeiChM5RWVWL/rIDYUlXqX5Uz8Eg9+bu6W1tv2HcbUBZvtKh4ABnkiorD97zfrMfyZBd6/Pf3k3170h6ntL399MZ6YvQ57SsvsKB4ABnkiorAt27o/ou1Ly8oB2DtiNipBXkSmichuEVmlWVZfROaIyEb3//WicSwiIqexswtmtGrybwIY4bdsIoB5Sqn2AOa5/yYiIi/7u1xGJcgrpRYA2Oe3eAyAt9yP3wJwfjSORUSUqBKvl7y9OfnGSqmdAOD+v5HeSiIyTkTyRCSvuLjYxuIQESWoRM/JR0IpNVUplauUys3Ozo53cYiIzIuw6h6LWSvtDPJFItIUANz/77bxWEREpMPOIP8FgCvdj68E8LmNxyIiirtwa+Z2Tm8WrS6UHwBYBKCjiBSKyLUAJgMYKiIbAQx1/01E5Fg/b95raf1YNNSmRWMnSqlLDJ46Oxr7JyJKRP5Bev/h43EpRzBxb3glInKKFE3U31VyFDOXFRquW15R6U3T2DniNSo1eSKiE5F/Dv5YRVW0vvS1X7Cl+BCGdW6CGtV8Q21ZeQU63v91LIrImjwRUbQ8N2+j9/HuA65Jxyp1qulHj1fGrEwM8kREcZYMc9cQEZ1wzNzuL953f2WQJyKyQaLMY8MgT0QUY/4Ntgk/nzwREemzM4CbwSBPRGQHC/mahJ/WgIjoRGRqrhqdCB7LfD2DPBGRDdjwSkREAABlY+KeQZ6IKExm0jV2DnQyg0GeElJewT6UV8Ru6DdROIINhpIgZ4Bgz0UbgzwlnOXb9uPClxdhypwN8S4KUVB689KEg7NQ0gml+KBrYqcNuw7GuSREwZm5SYh/AP8obxv++90mm0oUyBFB/vCxcmwsKkVOwxqok5Ue7+JQlMR7zg+iSBhlZO6ZsSKm5XBEumZjUSnGvPgTlv6xL95FISIbHGf7TNgcEeTJmRKlnzHF16rtJWg/6SvMXVMU76KEJd5XpAzylLDi/eNIJOt3HcTpj3yL3QeOxrsoMZe/bT8AYP763fEtiEWJUkmxPciLSIGIrBSRfBHJs/t4lPxi2Lssabzx0+/48/BxzFuXXIEumuI90ZednNC7ZpBSak+MjkXkWE4OdEaS/aRv52hWMxyVrvlpU+juTLGilMKUb9dj+/4j8S6Kj4/ztuH6t3lBlWw8gS7eoyfJPCsDnuw8kcUiyCsA34rIUhEZ5/+kiIwTkTwRySsuLg7rAJ436PWFv1varrJSYcGG4qBn2rLyCuwqsZ4HXV90EM/P34Sb3l1qeVs73T1jBeYkSQNWvGtAiSXJq7NR4dzvQ7LfNKSfUqo7gJEAbhaRAdonlVJTlVK5Sqnc7OzsGBSnyrSffscV05bg1R+3GK5z5/Tl6P3veaiotPYplFe41j9W4dwvJsVeop73pi38HVO+XR/xfhZv2YsZSwt9lpm5j2oii/dHZnuQV0rtcP+/G8BMAD3tPqZZ2/YdBgA8MXsdthSX6q7zzepdAMIfvpzcX09KFFXpmsT0yKw1eH5+5KM4/zb1F/zj4+W6zyXqCc6Ild++nWk4W4O8iNQQkVqexwCGAVhl5zHf+eWPkOtUVqqAVMAug65pkb71idxotKHoIJZt/TPexQgQz/csZ+KXePLrdfErgIEE/hrZLpF/Q8nA7pp8YwALRWQ5gCUAvlRKfW3nAR/4LPQ5pM19s001Ph4+Vu59bPV7Fu9aR4XOiczfsKcX4IL//hyjEpkXr/eu1xNzAQD//X5zfApghk1vzuodJfhgydaQ6932wTLc/9lKW8pgh+37j1hOtUZbvGOBrUFeKbVFKXWa+18XpdTjdhwnnJzd3LXB+xuv23UAnR/8JuwviOfyKyWMakjJ4eO4b+ZKHD1eEdaxAaDtfbNx24f5YW+fiP7Yewg5E7/Ewo3We+N+t243ftjg27CvlPL5fIsOlEVcRrtEkq5Ztb0En+dvD7rOqOcW4t5PQwfvL5bvwLu/hD4Z2MFqsNy69zD6TZ6P5+ZtxMrCEix3D6oKR8mR45i9cqelbaz89JO94TVhHQ8SwFdtP+Dzt153qD2lZXh94e+6NWbPIrMf9MrCEkyYno/KSoWn527A+4u3Yvqv2wLWm7+uCFv3Hja1z/9bvsPcwQ2UlpVj5LM/Ys2OA6FXDuHsKd/jvcWhU2nz1hYZBttfC1yppU+XFeo+H8zVb/6KK6ct8Vn2+Jdr0fa+2RHV9K56YwlyJn4Z9vZmRdL4OPr5hbjdfcLfsf8INhu0P8Xaul0H8POm0CfscF/5jhJX9+VFm/fi3BcWYsyLP+Huj5dj6R9/Wu65dceHy3DTe7/hj72HwiyNi6cdMJZO6CAfzG9BctXz1xVh277DmDA9H4/OWoM1OwODoNWwMe6dPMxcth27Dhz1fgH1vojXvJmHQVO+t7j30MorKlHpF+x+2bwXa3ceiEqvic3FhzBpZuhU2rVv5eG+ma4aZaSVm9/3BP9BvvlzAQBYDvKVlQrPzN2AvaVl+H59eN1+wxVpja/v5Pk4e8oP0SlMGA4cPY4NRa4ppEc88yMufW2x6W2tNk563yvNWeLjpYX460s/46UfrKXkCv90nTDKyq1PlKYt9znP/mh5+0id0EHeqIZQcvg43l/se0mqDbjXvJmHoU//gJIjxwFUdZfUW99yLt/EOhWVKuLa44t+81m3m/QVxrv79L+3+A/8eeiYpTLF0m4TaZX564ow6D/fY9YK36uZ2z5YhsdmrQFQdZV1vKISRRbmhPlly148M3cjJppIb0SLN10TQZSPd276f75eh1Mf/hbDnl5gabtwG16rUqaBz32cZ/1q0LrAAx8sK9dZz97fmGODvFIK+w8fC7rOfneQdm1Q9XD++tCDhY4eN3lGD/INnb1yJz5b5sqVxnoq1f/9JrB2PmdNEdbsOIBJM1fhzo/yoxJY9Ow+cDSsgON5Jxdu2oOn52wIuo+1O121xdV+qaYvlu/Aa+5Bc54U3N0zlqPXE/NMl6Pcfdwjx4zbTMorKqM69XU0Opjc9F58B+a9FGGDttWvoTdlqvPunUgD7RwR5PXi6Du//IFuj8zBluJS7D98LOT9Qo9XKhQfLMN1b+VhwnT9frpW6FwpBrjpvd9wx/R87Nh/BHtKXSckpRQW/x5ZcFi9oyTsbcvKXYFr3+HjYTX27Sw5EjS47Tt0DD2fmIcnZq8Nu4wA8Oy8jZizZpfh82UmGq09n803q62NADYzxcBz8zbiry8tCpr2C0ckocnK65y3tgi7DybGjJfhtkd47jC2aEvgdCdm3sfpv25F7yAn/3lrTb6fcT6fOOLOUHrmu2fr21B0EDe8+xsu6dkK//7LKb4rad58/0Y5K/Q+Q2+6xsT3s0CTO1YKWBfhbe82FlU1rP156Bjq1cgwva32tVj5cR0+Vo6hTy3wztWz+YlzcKy8ElkZqT7r/em+uvpu3W48MLqz6f3rKQ9Sk3/OxMCccNMAnvclWGVwrfsz9AQaK5RSKDlyHHWrV31unquOWFRAKysVrn0rD60b1sB3/zjL/gPa5I7p+YbPmXkf//mJu23IYOVr3wreDdta75okHQwVT573zNNT46O8wJ4qVhpyNhcHNuKZ+Qy163y2bLt+w4tmpb+9ssh0mYxoX9ehY/o5QMNtvZe45mzbdxh5BfuwdudBn8nY/vLSzzj5wciHRMxZU4Q9pa7P0H/6CU+wvffTlbhO84PzXI2EEm4NMcWbxgprcx8VlQo/b67qYbJj/xG0vnc2uj0yx9ZeMKG6VAKuhuvKSoWFG/foBqHyiko8P2+jz3iSaCvYc8hb6fl4aSGOhdHwqcfKbz/ZMzuODPLlFZXej/ChL1YDcP2YZiwtxOjnq4KslQ9v/DvWZm7U7nufuxHzjun5WLPzAA4fK/fpA7+3tKrtYEcYk6H5M6o9Wum+pa2FGL1PhX8explPfocLXw48MRn1SdY2Uh89XhE0jXb4WAWufzvPe5Xlf4XjKeMHS7ZirubS+fN8c11Hw63J73V/nto0gDZFNfSpHwwngVu9o8RnwrtXFmzGpa8u9vbh1/bFLtDpHeT/Uby+8HeMeWGh9++iA0fxzNwNUEohr2Cf4SyotwcZQ/HDxqoeQ28tKsBlry/G16sCU2Of/rYdU+ZswDNzNxruy2PfIeP2sVHPGfc4Oes/32PaT1UTDz725ZqQxzLDzG8/0qkkSo+WR7R9tDgyXdNu0lc4s33DgOX+c2JYuQFDwd7DKCuvwPPzzM3P4flgf9u6H90fnYP3r+/lfa7zg9/4rHvrB8t097GisATFB8vw3brduPiMlgHPHz1egcz01IDlpWWBNdmNRQcx1FSvBlfJl23dj39+ssLntfgb+Yy57mCbdlcF56fmuBp8t+w5hE4PfI2zOzXC61edobud5wSwde9hjHnxp4Dnl/y+D2/+VBD4CkyevVO8KRBz628sOojhzyzA0M6NA56b9lMBepxU37Xe7qoauOc8MmNpITo1qYXRz7sC8owb+uDClxehc9PaAKDbu6e0rBxKKYiIYSP4o7N8g97tHy7DL1v2YVDHRrjw5UVI1etaorF6RwlGPbfQZ1n+1v3ex3+4x2Ts9Kt8XP92nvc3dsigx4jWFdOMu0r6N44Hkx/BgCYtMx95iggqlMLOksATpZl2ryMWBjPaeSJwZJA3y+ql33u/bMUL34UO8qVl5d4fh8fKQuuNoZ8u245P3b1vujSvjS7N6vg83+mBr3H38I64tn9rb7AvPlim2+hoLsD72u2+IjAKgkbdwfyt0Lx2/x9psBPtb+5gc7CsXPfKwNPP3aPk8HGUV1aaSsP8/bVfUBqk/EopvL3oD4zp1sybG/94aSEqlfWGWiCwguEZqOY/xkI7Qvr2D/NRdOAo1u8qxdZ9wfv852/bj24t6+KIu9eXZ0K9UL2YvgvzTlNz1hShvvt9MTrE0eMV2LS7FF2b18GGXb6pJ72TtlbOxC9RMHlUWGWLFs8nMfSpBWhRL8vnOf8TYyJzZLoGsCePtmCj/sAX/yD491cDZ9KLdJIloy6b//vNetyvma/njMfn4pUFVbnrC19aZGl6BLPvW7DLb3+fadIndk4dcNoj36LHY3NNrRvqBjP52/bjoS9W454ZK7zLwv0Ipy4I7Dr41iLf0b+efft/T56YvQ6f/FboHe1r5PwQQVPPaz9uCaihA/CZw8aT4tP7WmzZ4wrcRpWATg98jdHPL3Slp/xel/9JO1h3VD17S8vwlt9JfktxKT79zVz/d6VCz+3k7SobwfQirmNFtHnEnBvkbbgA8h/duNxdQ91SfAibNJfoy3Vq7ZHOiR2sH/2MpYX40eAEtOvAUSzTXH7rCfUD06sN3m0wHayeFYXBjx9t93xSFZitpAK0PCejRZvN3W2sqjumb+5aRPDE7NCzWnrSNaHmOormSfixL9fiPb9Bf0ePV3iv4ICqK63XdO654DnxaNNTCzfuwejnf8TBo1VjUErLykN++298bykK/zTXZnSsvBIXvbwID32xGjM1U1wMnvID7vxouamG4B0lRzH8Gf2r210lRy2fdIIxE4vsHCfj2CC/NYZzRNz18XIMeeoHXDltCeav07+UD5YaMGPyV+sCRqlqXf76Ep9RqlqhbkH4yKzV3sfak5XHz36Bbu3OAwFpFqPYdMv7v4UMTEePV9jWhWzfofCuHDxXbQfLyqt6Zpk4T49/J7wBR//5doPrEFGaVjdU9z6rdpYcxd5S/fdSWwm47PXFWLX9AE55+FvvMpHQr+v79cXo/z/f+SwbrpNiXFFYgrs+Xo4t7kZpvTEtCzfuMdV1dUNR1ZXIul0HvL2yev97XkDPsI06v4to0mtbihbH5uS37Yv9vVV/2FAcMNOhx7PzQvdACCZ/2/6QjU5G82p8FWL2vA+WVHUvfddgErGdJUdQqYCGNTMwUqcbqFHj26wVO5GZHrwu0ekB+2af1ptywl+oNe6ZsQL3zFih25gfyvpd5q8knp+3EXWrpwddx+57vL69qMDwOaNUWKjvZbhz5awv0h8vEmrivXEWT7TnvfATVm53XX1fnNvC0rZmeOovv+gMyvI4FMUrB3+OqMlHmjNzit7/1h+dZ6UXkf/smx4bikrRb/J83PWRfprm8teNB5OZngLCBut2HQz64wL0UyD+cxcBwI9Bpjg2unOYp4ZuxpQ5G0KGcO1h/FMKZscHBGMmteQ0ngAPAB/ZNKdNWXkFxk79xZZ9h+KIIG+lEZDC45mhctYKa3NqJ4JY/Lhmr9yFcSZuRBPKNAs3o/c/qXe839b78RDCa0RVgO7Vr5bZ9ohwOCLI8+5g9rv6zV/jXYSE963BACgrCkzeK+Do8QrvLKg+4t2Vw+EmTM+3PANsZaXCFp0R81rLtu4P2uYWCUfk5MO5+xJRspowPd9whtVwexMls+/W78avYU7q95TFeyXo3TsilDOf/C70SnB1h755UDvL+w/FEUGeMZ5OFP/+KnjOPNikbU519RvhX2Wamcgu2dmerhGRESKyXkQ2ichEO47BmjwRkT5bg7yIpAJ4EcBIAJ0BXCIikc0vq3ucaO+RiMgZ7K7J9wSwSSm1RSl1DMCHAMZE+yCsyRMR6bM7yDcHoJ3IvdC9zEtExolInojkFReHd1PkLQly93kiokRjd5DXq2L7tAwppaYqpXKVUrnZ2dlhHeRolG4kQETkNHYH+UIA2onQWwAwd0cHIiKKmN1B/lcA7UWktYhkABgL4Aubj0lERG629pNXSpWLyC0AvgGQCmCaUmp1iM0sY7MrEZE+2wdDKaVmA5ht5zHaN65p5+6JKIGt/tdw/LH3MM4Jcq/YE5lD5q5hXZ7oRFWjWho6N6sd72IYmnfXQFPrzbq1vy3Hd0SQt3uObaJwjR/YJqr7m3LRaZh310D83y32BIQTzdpHRlje5s6hHSyt37xuVsh1LunZCl2b1wm5XjgcEeSpSs1q9mbgamda33+/dg1sKEly8NzsOloUgLbZNXFKizpIT/W9gu3avDb6tInve/3M37pZDoLR0ja7huVtsjJSff5u0zD0Pm4dHP1JxFJsTEYwyMfI5L+cEpPjrPrXcFv3P+vWMy1v06Ju9ZDr1MlKR8HkUT7LsmtVs3wsPTNu6IO7h3fEXREEn6X3Dwlru14Wg+6oU5uaXreZXw3xb2e0Qnqa9Z+0//tulfbEcv7pzXHb2e0j2l+4JMyR7wv/Ocj7+JMb+9pynOFdGgd93s5R+44I8rGYQvsWnSlAcxoEBq/bDb7gY3u2inqZ4qGVzmsOpXGdzJDrjOnWLGDZkJOD/zDMys2pj5sHtcPZFvb3j2FVJ4TxA9qgQc3gJ5y+bRvofvZ1s9Itza3UvVU9AMDADqEHBk4f1wcvXtodax8ZgdvPbo+xZ7QMuY2H52RSv0bwK41GJk60enfF0t7G8OMb+pguVyTCDZMN3Z/tjWe1Rb0aGfjXeV3w4bje3uev6dfa+7hDmJ08nr+kO357YKjh833a2ncF5oggHwt6Nax2jWr5/P3pTX0xIcqXqmd1DG8UcDSd7xeAR1uobQK+P74W9QLzk+9e2wuTRp0MAFgy6Wz8pbtn5ovgZ+92jQJ/cJf1boVW9fVPRFYa505q4LpsH3VqU9x7jqtszYKcrN6/vrfhZzX3TlfDm5lUl+e9MlNvaVInE6NObYqsjFRMGNoB6akp3huiX9nnpKDbPjamK4DAAL3p8ZHexyO6NMG71/Xy/p1qkFPQq2TlPzgM3VvVBRA6+H55m377gtEVhtF9WMOtDGemp2Lj4yNxz/COAIAr++agt+bqRPu5/uei08I6RkZaStATav8w7h9sliOCfCxq8tVCXAaf2b6htxYWTdG4jJt2VW5E2/u/vS9c2j3gB/jo+V19/n52bDfv4+vObI2Lc1tg5cPDMPOmfgBcNwT36N++IaqluXKjjWpl4rQWdQEAlSFmq2iXHRjkz2yfjfc0gSkcGan6n/XgkxsF3c7oMr5tdk18O2EA5v/jrJDHbuo+kRg11imTX/ZQVy11stJx7mnN8OoVvt+NNM1rr5mZ5q3lAoFXrie5/zZKq3lKKoKgn0mXZnUw986BeOXyHkHLDLjy7k9eWBVotTXrFy/tjotzW6BRrWqm0l5zJgzwPk5PTQn4/DwntYy0FAzu5Prsa9jc5mUHRwR5M05tEVnLtX6wdX2Np17eA+9cG1lg0XNqizo+DTKhGmeu7e+6rOziV2MN1cX07Wt6hlU+rYHtfWuxtbOqLtdrZabjyQtPQ63MdFM9oc5012o8KRyjoHtmh4Y++eAHR3fG0JMbG95UW+v1K41PfGmpgt5tGiBFfC/V61loRPUEQE9JOjSuhYY1q+HekZ3w/vW9vEHohoFtfbYb0bUJ3rj6DFzTL8f0sbQ8PTQa1TZOs/w8cTBSUgTPX3I6zsipH/C8J6BNGNoB9WtkeL8f2iB4equ6+HBcb6SmiOHdjKo+BkG/dsFrqu0a1cTwLk28f8++Tb/tx1OGxy9wVSreuLrqu9u+cS08eeFpWDJpCF68tHvQ43VoXBPtG9cKuo72SuTpi7vh2bHd0FanYpHoTpggP6yztfzuZb19c+jBKtT+NQBPML5nREfd9c22zo/o2sS77wt7tMD8u87yeb5g8igUTB6FWbf2x4OjO+OB0a6p+v1jnDY/ql/+wGXN62bhyQtPDbqdtjbvH7w9u+ypE0R81wjUJrsmCiaPQp+2DdCpSS089Tf9S+QaGWn4YFxvfHxDH8y8qS+u6d8aKSmCWplVr1ebWwWAJ/96Kl67Ile3puv5XJRy1U63/HsUepxUdXV2y+B2+Nd5XXB2J/0a/SkmusCNH9gWfds29L56bU30/lEnQ0QwqGMj7+ce6rPzd9fQDph1a390alIbz47thneuDTyBh7o6nHbVGSiYPMp7NdG4tuvqQrtVs7pZaFonC5ufOMcwDZbm/iGE03PEaJ+eff2910k+ZdQTKmUVysPndUHP1vVxWsu6qFM9HWO6NQ+9UQJyRJA3cwEb6qz90LmdcYmmcfSeEZ18ntf7YRhVGD212OvP9O0j/eG43riqbw4u7dUKdbL0f7wZmrRQioj3S312p0bIaVgD947sFLBN1+Z1cE3/qhqntliX9GyJ01vVC5pH76h5bzz58OoZqQE16NN0roau6dcaN53VNmC5iGDdoyPw/vV+VzgWUmsigq/vGIDRpwY2ymqdkVMfp2tSZdrcZ2+/3i0Xn9ESQ/xO+J6uiJ5au9HVRrW0VFzZNweT/3oqxg0I7P/uubyv7tctT48niKemCNY+MgJbnjgH1/l8X5T3tVzVN8dvqbG01BRvbX5Mt+Y4s31gO4HVcSWe9bU/ATNx+5mx3XBd/9be9Fs4avi9l1YGPtYLkgM3k/Xq0qwOPhrfB5npgZ+nXttSonJEkPf39R2+l3ofjuvtcymop1+7hvDEtEfHdEHtzHSfPHywXJz/127mTf3w8LmdkZ6agntHdkKnJq4g2rtNAzx8Xhc0rZOF5Q8N092XNuAK4K2VenKlRicHLW3edoD7R/7s2NMN129Uu6pB0ZPySU0R9GrjqoX/vddJ2PT4SHzqzqdrPXhu54AToqecmempPjlewDdPmyj8A0eoAJBdqxruczfG+pt750AsuGeQ7nO+x6ySlZGKFIPqrsBVo7ywh35jYyx43g/t+zR+QOCJ3V+LetVx/+jOhq/NjF/vH4K7h1ddEec0NN+7y86R8J/f3C9pBqQ5Isj7N0Z1alIbMzTdtvxrc+9e2wsL7vb9IQqA1g1dl85N6rjO0p7L70fHdAnZ1UyrdcMauMpdKxw/sC2+vmOA7nqL7zvbp/HHn4griP5jWAdviuCiXPPd5ICqoJqaIkH7nf86aYhPFy8RQdM6WSiYPAo9W9dHWmqKYe8KwDcwDjm5kW6tXysaPz8zNWYzPDV5z4nHbD33o/GBXQPbNarp01gZitEJpXXDmhh9atOgJ+dwWe2o4Lky7dayrnfZKRG2cRmZenkPn/EM1TPSvAP8mtfNwpSLu5nel50ViQY1q9n2HkSbI4K8nlzDXLCrN4e2v3eKuH6cV/fNwfvX9cJQv8v5Tk3184ORduppXDsT7RvXQv6DQ72NXdrAJRDUzkzHLYPbe2tDqSmuNMi6R42HYxs1PGoHffjLrlUN9WtkhN1TybPZSQ2q47UrzzDsaRLu/v0Hkz12fteAz0nrjavP8KkBBjPjxr64Y0j7qstyk2Xs2dr4OxaKp3HZaHK91BTBC5d2t2Wou5mGaa3mdbMw69b+eOT8LlErwx1D9MeTDOvSBLcajDUZ3KmRpRHdGWEMDIsmveN/elPowVbRlnz9gSz4aHwf7C0t8/798Q19sHp7ScB6E4Z0gIhABOir6QUQ6qfguYKItMZQt3oGLunZCvPX7fbJ//kPufbQyxH6lkt/uaebohlWX5L3vTC7f/eKtUxOkzC2ZytM/HSl9+/LegdvVBvUsREGdQze5dGjY+NaOLlpbZSVVwCIzVxIF+W2xJCTGwfNG2t5GhgbWLiijKZon2zuGNIBz8zdaGrdcH9fV/XNwd7SMrz64+8Bz43sGjx9axc7ulmH4ugg71/TOiOnvm6XsUgHIkTjstCTCclMT8G8uwbii/wd+JuFEYxauTn1sXF3fO57G2rItzaAfnJjH7SoZz7H2qJeFgr/POJt44iUiO8JMS3FVfPq1Try0YdmvhJmAzzg6vnTuVlt7xVfOJrWycTOkqMxGVdiB6sn38z0VEwa1dkb5K/r3xqvLfwdaSmCO4bEZ36deHBEuibS76zeyMlYG9ghG9f2b43HLzgFbbNrekcwhqNmtaoau9EP+mqDftjeRjab8pnaRrweJ9X3ds8zwzPasLaJxudwpKYIvp0wwNSgHI+5dw7EB9f3Dlh+62BXyqFxkP7qVqSlpmB4lyYhT6J6Ftw9CB+O623r/CjhMDM7IxC9mwJ5ei/VyUqPqDE42TgjyEcY5Wtk6F/Q+KcgZtzQB//U9CSJZoUoLTUFD4zubKnRzoj20tq/9uOpCT50bvD8qtV4YPW9SLB449WhcS1LoxrbNaqpO+/IX3u0QMHkUahu8N2KpVYNqgd0PoiH6X5jFubdNRCrLUyoF+nv3NNxIFG/e3ZxRJDXalnffP/V804L3v/aX25Ofdyo1yc8wW5aEmxir1evyMX6x4wbbT01T6uTg1XV0EOsZ2mvsZGIZbKDt/dQnF6w/4ycmemp5k6o7oJHUux61c2Nto6mRIkK8a9mRNnnN5vvuzrl4tPw0LnG/XhD9elO1Nymtlj+ZUxNEaSmGDfANqqdifwHh6J2ptWUiMkoD0ur+ejeqh5GdGmCf5jsNWOmDAn6EdoiWWuwkRZ72QNDkZ6WgooK16d9skFvOadyXE3eSn/29NSUkFPIuuh/zTy54Xh31QrXoI7ZujnjutUzws5ZhqzJR3BmzEhLwcuX90iINpRkdHVf19iN+jUDfyN3Du1geQqFWAv3q1OvRgZqVktDnerpmD6uN/779+Dz2jiNbTV5EXkYwPUAit2L7nPf1NsG9tTHQn2pHju/K05pXht9bZwL2k7ayZ0i1bJ+dbTJroGHzwue669q2E3SamUSu6Z/a5/pL7RuO7u9qZt9/HjPIKSlhv/ZTbsqF3PWFFnaJppfFas3cYmGeF/x252ueVop9R+bj+FlV7A1+pLVyUrHOBPDu2NNO99HrL5f1dJSAyZQI+dpaTBXv1mDOzXG4E7h3gwmuZJrVSOo41vu5Mwz+PGcKc0OrDG936juLXZExPKNPU5EL1/WAz1z6tt6f02KXGqSXvUlSocMu2vyt4jIFQDyANyllPrTfwURGQdgHAC0ahXeLfKqu1vo69eITp9kf4nxUYUnkhy4nRLhdzusSxMMCzFxHcXf+ac3x4rtJbh7WHQa3GMhkSoOEdXkRWSuiKzS+TcGwEsA2gLoBmAngCl6+1BKTVVK5SqlcrOzw7vV3YD2DfH4BV3xwGj9mQHD1dU9p3VdnZtFWOmqGQ+eq5pQd7SKtQQ951ACy0xPxRMXnGJphHCsNDVx/+J4i6gmr5QydQt7EXkVwKxIjhVi//h7r8huEKDnwXM748IeLdC6YQ2f5d9OGIDsKAxastN955yMFvWqY1jnxKqp6s1NTslv3IA22P7nkXgXIyHFu2JjZ++apkqpne4/LwCwyq5j2aVaWqrPzSg8OoS4AUkiqJWZbnhbtkSQKPlKig6j+fVPZIlSkbEzJ/+kiHSDq/2yAMB4G49FScIzs2a0JhkjouBsC/JKqcvt2jclr0a1MvHR+D4BNxsnSkYNamZgZ8nReBcjKMdNa0CJL5KbbRAlkgY29eiLpsTqekFElETM5N3j3aGMQZ6IyAZW2l3tbKNlkCciCpNRcI537V2LQZ6IKEzJMNEegzwRkYMxyBMRhclMPT7e80cxyBMR2SBRUjkM8kREYUqQOB4UgzwRUZS8eKnr1oLxnpRMi0GeiChsvlX5EV0Ta9ZXgEGeiMhW8a7UM8gTEYUpWE4+UdL1DPJERA7GIE9EFKZ/ndcFzesm9q1AGeSJiMLUrG4Wnrr4tKDr+Pe0ubZ/axtLFIhBnojIDgZJ+QdGd8ayB4bGrBgM8kREEWivuedzojS2ajHIExFFoH6NjHgXISgGeSIiB4soyIvIRSKyWkQqRSTX77l7RWSTiKwXkeGRFZOIKEnFeTRUpDfyXgXgLwBe0S4Ukc4AxgLoAqAZgLki0kEpVRHh8YiIkkKi5OcjqskrpdYqpdbrPDUGwIdKqTKl1O8ANgHoGcmxiIgSWZ2sdNPrxnL2ykhr8kaaA/hF83ehe1kAERkHYBwAtGrVyqbiEBHZZ8Hdg1ArMy0hpx4OGeRFZC4AvanVJimlPjfaTGeZbmZKKTUVwFQAyM3NjfdcPkRElrVqUB1A/O8CpSdkkFdKDQljv4UAWmr+bgFgRxj7ISJKairOLa92daH8AsBYEakmIq0BtAewxKZjERElHEfc/k9ELhCRQgB9AHwpIt8AgFJqNYCPAKwB8DWAm9mzhogo9iJqeFVKzQQw0+C5xwE8Hsn+iYicSGLYwZIjXomIHIxBnogoSrR5+AlD2gMAqmfY1VPdnPgenYjIoa7q1xpX9Yvt3PF6WJMnInIwBnkioliLYe9KBnkiIgdjkCciijM7B04xyBMRORiDPBGRgzHIExE5GIM8EVGMxXLuMgZ5IiIHY5AnInIwBnkioigb3KlRvIvgxblriIii6OeJg1G/Rka8i+HFIE9EFEXN6mbFuwg+mK4hIoqxWN4YkEGeiMjBGOSJiBws0ht5XyQiq0WkUkRyNctzROSIiOS7/70ceVGJiMiqSBteVwH4C4BXdJ7brJTqFuH+iYgoAhEFeaXUWsDeaTKJiJwmljHTzpx8axFZJiI/iMiZNh6HiIgMhKzJi8hcAE10npqklPrcYLOdAFoppfaKSA8An4lIF6XUAZ39jwMwDgBatWplvuRERBRSyCCvlBpidadKqTIAZe7HS0VkM4AOAPJ01p0KYCoA5ObmKqvHIiIiY7aka0QkW0RS3Y/bAGgPYIsdxyIiImORdqG8QEQKAfQB8KWIfON+agCAFSKyHMAMADcopfZFVlQiIrIq0t41MwHM1Fn+CYBPItk3EZFTcVoDIiKKCgZ5IiIHY5AnInIwBnkiohjjjbyJiBwsKz0VNwxsG5NjMcgTEcWYiGDiyE6okZFq+7EY5ImIHIxBnojIwRjkiYgcjEGeiChOsmKQk4/0zlBERBSmj8b3wTeri1Czmn2hmEGeiChO2mTXxI1n1bT1GEzXEBE5GIM8EZGDMcgTETkYgzwRkYMxyBMRORiDPBGRgzHIExE5GIM8EZGDiVIq3mXwEpFiAH9EsIuGAPZEqTixlsxlB1j+eGP54ycRyn6SUipb74mECvKREpE8pVRuvMsRjmQuO8DyxxvLHz+JXnama4iIHIxBnojIwZwW5KfGuwARSOayAyx/vLH88ZPQZXdUTp6IiHw5rSZPREQaDPJERA7miCAvIiNEZL2IbBKRiXEuyzQR2S0iqzTL6ovIHBHZ6P6/nua5e93lXi8iwzXLe4jISvdzz4mIuJdXE5Hp7uWLRSQnimVvKSLfichaEVktIrcnWfkzRWSJiCx3l/9fyVR+9/5TRWSZiMxKwrIXuI+bLyJ5SVj+uiIyQ0TWuX8DfZKp/IaUUkn9D0AqgM0A2gDIALAcQOc4lmcAgO4AVmmWPQlgovvxRAD/437c2V3eagBau19Hqvu5JQD6ABAAXwEY6V5+E4CX3Y/HApgexbI3BdDd/bgWgA3uMiZL+QVATffjdACLAfROlvK793kngPcBzEqm7457nwUAGvotS6byvwXgOvfjDAB1k6n8hq8rFgex9QW43sxvNH/fC+DeOJcpB75Bfj2Apu7HTQGs1ysrgG/cr6cpgHWa5ZcAeEW7jvtxGlwj7cSm1/E5gKHJWH4A1QH8BqBXspQfQAsA8wAMRlWQT4qyu/dZgMAgnxTlB1AbwO/++0uW8gf754R0TXMA2zR/F7qXJZLGSqmdAOD+v5F7uVHZm7sf+y/32UYpVQ6gBECDaBfYfSl5Oly14aQpvzvdkQ9gN4A5SqlkKv8zAO4BUKlZlixlBwAF4FsRWSoi45Ks/G0AFAN4w50ue01EaiRR+Q05IciLzrJk6RdqVPZgr8n21ysiNQF8AuAOpdSBYKsalCVu5VdKVSilusFVK+4pIl2DrJ4w5ReR0QB2K6WWmt3EoBzx/O70U0p1BzASwM0iMiDIuolW/jS40qwvKaVOB3AIrvSMkUQrvyEnBPlCAC01f7cAsCNOZTFSJCJNAcD9/273cqOyF7of+y/32UZE0gDUAbAvWgUVkXS4Avx7SqlPk638Hkqp/QC+BzAiScrfD8B5IlIA4EMAg0Xk3SQpOwBAKbXD/f9uADMB9Eyi8hcCKHRf+QHADLiCfrKU35ATgvyvANqLSGsRyYCrQeOLOJfJ3xcArnQ/vhKuXLdn+Vh3q3trAO0BLHFfFh4Ukd7ulvkr/Lbx7OtCAPOVO8kXKfexXgewVin1VBKWP1tE6rofZwEYAmBdMpRfKXWvUqqFUioHru/wfKXUZclQdgAQkRoiUsvzGMAwAKuSpfxKqV0AtolIR/eiswGsSZbyB2V30j8W/wCcA1dPkM0AJsW5LB8A2AngOFxn7mvhyrvNA7DR/X99zfqT3OVeD3crvHt5Llw/ks0AXkDV6ORMAB8D2ARXK36bKJa9P1yXjysA5Lv/nZNE5T8VwDJ3+VcBeNC9PCnKrzn2WahqeE2KssOV017u/rfa8ztMlvK7998NQJ77+/MZgHrJVH6jf5zWgIjIwZyQriEiIgMM8kREDsYgT0TkYAzyREQOxiBPRORgDPJERA7GIE9E5GD/DxgyVNK6SsftAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(res.fun)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see much better picture of residuals now, with the mean being very close to zero. There are some spikes left. It can be explained by outliers in the data, or, possibly, the algorithm found a local minimum (very good one though) or didn't converged enough. Note that the algorithm worked with Jacobian finite difference aproximate, which can potentially block the progress near the minimum because of insufficient accuracy (but again, computing exact Jacobian for this problem is quite difficult)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "## 1. SfM pipeline \n",
    "\n",
    "To get the context of below questions, take a look at the code above: The same questions have been asked at different places above as comments in the code.\n",
    "\n",
    "1.  **Basics** - How do we know this (`camera_ind`) information in practical setting? In other words, how do we know observations in `points_2d` belong to which camera. Explain. \n",
    "    - Ans-1 - Basics:\n",
    "        * We know the ordering of camera views and are able to extract features from the same. \n",
    "        * This ordering of camera views along with the ordering of poses (which are also known) allows us to generate the array 'camera_ind'\n",
    "\n",
    "        \n",
    "        \n",
    "2.  **Basics** - How do we know this (`point_ind`) information in practical setting?  In other words, how do we know observations in `points_2d` belong to which 3D point. Explain.\n",
    "    - Ans-2 - Basics:\n",
    "        * We use feature extraction algorithms such as SIFT (Scale-Invarient-Feature-Transform) and SURF (Speeded-Up Robust Feature) extraction to find noteable features such as a sharp edge in an image.\n",
    "        * We look to find the corresponding feature in another image plane. We then end up with a stereopair which has 2d-2d correspondances. \n",
    "        * The depth is estimated via triangulation, as the extrinsics as well as the intrinsics of the camera's are known, we can accurately recover the 3d world point.\n",
    "\n",
    "   \n",
    "   \n",
    "3.  **Transformations** - `project()` function: In the `project()` function, would it make any difference if I do translate first, then rotate? Why/why not?\n",
    "    - Ans-3 - Transformations:\n",
    "    Transformations: If we do a translation first, our transformation will look like: $R(X + t) = RX + Rt$. Hence our translation vector is no longer independent and is coupled with our rotation matrix. This is in violation of the assumptions made in the other parts of the code, such as the Jacobian formulation which assumes a rotation first model.\n",
    "        * Yes, It will make a difference. The order of these operations is not permutable.\n",
    "        * Our translation and rotations are described according to a certain coordinate frame. If the translation is done first then the rotation will yeild a different result to the one we would get if we rotated first and then translated.\n",
    "        * In order to obrain the same results by translating first and then rotating, the rotation center has to be translated by the same translation vector.\n",
    "        \n",
    "        \n",
    "4.  **Jacobian** - `bundle_adjustment_sparsity()` function: m above is not \"M*N\" (*2) unlike our lecture notes. Why is that so?\n",
    "    - Ans-4 - Jacobian:\n",
    "        * In the in-class discussions and lecture notes, m denotes the number of poses, that is, the number of views. 'n' denotes the number of points obtained from each view.\n",
    "        * In accordance to the above convention, we would get m*n*2(one for x, one for y) equations. This is the total number of observed 2d points captured in the process.\n",
    "        * Here, the points obtained from each image isn't fixed like we had in class (or in the notion page). We are given total number of points rather than the points per pose.\n",
    "5.  **Jacobian & Parameters** - `bundle_adjustment_sparsity()` function: \n",
    "    1.  Why are we doing `n_cameras * 9` here instead of `n_cameras * 12`? Recollect: Every individual motion Jacobian was (1*)12 in our lecture notes. \n",
    "        - Answer to 5A - Jacobian & Parameters:\n",
    "            * The reason for multiplying by 9 here and not 12 as discussed in class stems from what that number signifies.\n",
    "            * The reason for having n_cameras * 12 in the lecture notes is that we are solving for all 12 values of the projection matrix P which has dimensions 3x4.\n",
    "            * The reduction in parameters is down to working with the axis angle rotation convention which has just 6 parameters. The number of camera parameters is 3, one for the focal length and two for the distortion parameters along x and y.\n",
    "    2. Ignoring the scale parameters, what was the number of unknown parameters in our lecture notes in terms of `n_cameras` and `n_points`? What is it here in the code? Is it different? If so, what is and why? [Link of notes](https://www.notion.so/Stereo-Structure-from-Motion-9fdd81e4194f4803ac9ba7552df56470).\n",
    "        - Ans 5.2 - Jacobian & Parameters:\n",
    "            * Ignoring the scale factors, the number of unknown parameters solved for in class is as mentioned above where n_points is computed as the number of poses times the number of points observed in each view divided by the number of times each point is seen.\n",
    "            * The number of unknowns in the above code is (n_cameras * 9 + n_points * 3). The 9 comes from the new number of parameters to solve for in every pose. n_cameras and n_points are known a priori thus resulting in a different number of solvable parameters.\n",
    "\n",
    "6.  **Sparsity, Residual Vector & Jacobian** - `bundle_adjustment_sparsity()` function: Explain what you understand from above 6 lines of code by coding a simple toy example yourself to illustrate how it is different from what you've learnt in class. ([Coding toy example + elaborating in words]- both are compulsory.)\n",
    "    - Ans 6 - Sparsity, Residual Vector & Jacobian:\n",
    "        - The above code not only is a new implementation but changes the structure of the Jacobian itself. On analysis of the above code blocks, we can say:\n",
    "        - The structure of the jacobian is observed to be sparse. The above code aims to populate a spare matrix which dictates the sparsity of the jacobian. The above implementation sets the non-zero entries to 1. This sparsity in the structure of the jacobian is something that is wanted to be exploited when performing computations. The scipy function leverages this fact for us and speeds up the computation.\n",
    "        - If the number of unknown parameters is given by 'n' and 'm' gives the total number of points, the size of the jacobian as well as the size of the residual can be computed. \n",
    "        - The jacobian is of size (2m, n) it has non zero entries only for the parameters that are used in the computation of a specific point in a specific image.\n",
    "        - The size of the residual is: 2m\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "       [1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "       [0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "       [0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "       [0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "       [0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "       [0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "       [0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],\n",
       "       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],\n",
       "       [0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.]])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "camera_ind_exp = 10\n",
    "m = camera_ind_exp * 2\n",
    "B = np.zeros([m,16])\n",
    "i = np.arange(camera_ind_exp)\n",
    "\n",
    "practical = True\n",
    "# theoretical setting\n",
    "if practical == False:\n",
    "    B[2*i, np.array([0,0,0,0,0,  1,1,1,1,1])]=1\n",
    "    B[2*i + 1, np.array([0,0,0,0,0,  1,1,1,1,1])]=1\n",
    "# practical setting\n",
    "else:\n",
    "    B[2*i, np.array([0,2,5,6,2,  1,9,14,7,11])]=1\n",
    "    B[2*i + 1, np.array([0,2,5,6,2,  1,9,14,7,11])]=1\n",
    "\n",
    "\n",
    "B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Initializing R,t and 3D points for SfM given 2 images (`4 mark`)\n",
    "\n",
    "Using OpenCV functions, mention how you would initialize R,t (poses) and 3D points for SfM given 2 images and K matrix. You don't need to implement it, just mention function names with input/output arguments clearly and briefly explain what they do (You don't need to give detailed answers). A sample answer could be as follows:\n",
    "\n",
    "**Ans 2:**\n",
    "\n",
    "1. **Feature mapping** <br>\n",
    "\n",
    "* We use `cv2.SIFT_create()` to create a SIFT based feature detector.\n",
    "*  We then get SIFT features which are keypoints and descriptors using `kp, desc = sift.detectAndCompute(img, mask=None)` for both images as (kp1, desc1) and (kp2, desc2)\n",
    "\n",
    "2.  **Find correspondences from features** <br>\n",
    "To find correspondences from features we find matching features using a knn based feature matcher. The steps for it are as follows: \n",
    "\n",
    "* We use `cv2.FlannBasedMatcher()` to create a FLANN, with initializion of search =  50 checks and kdtree = 5 trees.\n",
    "* We then use `matches = flann.knnMatch(desc1, desc2, k=2)` to apply a flann matcher with k=2(since correspondences between 2 images).\n",
    "* Now we can get good matches of the 2 nearest neighbors which are p1 and p2 using ratio test.\n",
    "\n",
    "3. **Find fundamental matrix** <br> \n",
    "Now, we find fundamental matrix $F$ using \n",
    "`F, mask = cv2.findFundamentalMat(p1,p2,cv2.FM_LMEDS)`.\n",
    "\n",
    "4. **Get essential matrix E**  <br>\n",
    "Essential matrix can be calculated as follows since F and k are known.\n",
    " $F=K^{-T}EK^{-1}$\n",
    " \n",
    "5. **Decompose E to R, t** <br>\n",
    " Essential matrix $E$ can be decomposed into $R$, $t$ with:\n",
    "    - `retval, R, t, mask = cv2.recoverPose(E, p1, p2, K)`\n",
    "    where  $E$ is Essential matrix, $K$ matrix and p1, p2 are the corresponding points in two images.\n",
    "    \n",
    "6. **Use triangulation to initialize the 3d points for SfM** <br>\n",
    "Now, we use triangulation to find the position of 3d points and initialise them for the LM/GN in bundle adjustment problem/Sfm. (since the camera matrix $K$, point correspondences in the two images (p1 and p2), relative position of camera centers using $R$ and $t$ were known)\n"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}