diff --git a/multiview_structure_from_motion.ipynb b/multiview_structure_from_motion.ipynb index 8f5d649..e56f118 100644 --- a/multiview_structure_from_motion.ipynb +++ b/multiview_structure_from_motion.ipynb @@ -33,10 +33,184 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 36, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import cv2\n", + "from structure_from_motion import compute_sift_and_match, u_to_x, Camera_2, Camera_Rodrigues, triangulate" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": {}, + "outputs": [], + "source": [ + "files = ['pictures/DSC03919.JPG', 'pictures/DSC03920.JPG','pictures/DSC03921.JPG']\n", + "i = 0\n", + "f1 = files[i] \n", + "f2 = files[i+1] \n", + "f3 = files[i+2]" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "metadata": {}, + "outputs": [], + "source": [ + "i1 = cv2.imread(f1)\n", + "i2 = cv2.imread(f2)\n", + "i3 = cv2.imread(f3)" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "metadata": {}, + "outputs": [], + "source": [ + "u1, u2 = compute_sift_and_match(i1, i2)\n", + "x1, k_cam1 = u_to_x(u1, i1, f1)\n", + "x2, k_cam2 = u_to_x(u2, i2, f2)\n", + "E, inliers = cv2.findEssentialMat(x1[:,:2],x2[:,:2],np.eye(3),method=cv2.RANSAC,threshold=1e-3)\n", + "inliers = inliers.ravel().astype(bool)\n", + "n_in, R, t, _ = cv2.recoverPose(E, x1[inliers,:2], x2[inliers,:2])\n", + "P_1 = np.array([[1,0,0,0],\n", + " [0,1,0,0],\n", + " [0,0,1,0]])\n", + "P_2 = np.hstack((R,t))\n", + "P_1c = k_cam1 @ P_1\n", + "P_2c = k_cam2 @ P_2\n", + "real_world_coords = []\n", + "for xx1, xx2 in zip(u1[inliers, :2], u2[inliers, :2]):\n", + " a = triangulate(P_1c, P_2c, xx1, xx2) # x, y, z, w\n", + " real_world_coords.append(a[:3] / a[3])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now that we've computed real-world coordinates from image 1 and 2, we can move on to trying to fit a projection matrix for image 3 that minimizes the misfit between real world coordinates from image 1 and 2 and real world coordinates between image 2 and 3.\n", + "\n", + "\n", + "First we find the projection matrix for image 3 assuming image 2 is at the origin. " + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "metadata": {}, + "outputs": [], + "source": [ + "u2_p, u3 = compute_sift_and_match(i2, i3) # u2 prime and u3\n", + "x2_p, k_cam2_p = u_to_x(u2_p, i2, f2)\n", + "x3, k_cam3 = u_to_x(u3, i3, f3)\n", + "E, inliers2 = cv2.findEssentialMat(x2_p[:,:2],x3[:,:2],np.eye(3), method=cv2.RANSAC,threshold=1e-3)\n", + "inliers2 = inliers2.ravel().astype(bool)\n", + "n_in, R, t, _ = cv2.recoverPose(E, x2_p[inliers2,:2], x3[inliers2,:2])\n", + "P_2_p = np.array([[1,0,0,0],\n", + " [0,1,0,0],\n", + " [0,0,1,0]])\n", + "P_3 = np.hstack((R, t))\n", + "P_2_pc = k_cam2_p @ P_2_p\n", + "P_3c = k_cam3 @ P_3\n", + "append_row = np.ones((4))\n", + "P_3c = np.vstack((P_3c, append_row))\n", + "initial_pose_guess = P_2c @ P_3c.T" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We then multiply the projection matrix that we computed assuming i1 as the origin with the projection matrix we computed assuming i2 as the origin. This gives us a decent initial guess for the position of camera 3 relative to camera 1." + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": {}, + "outputs": [], + "source": [ + "X_gcp = []\n", + "u_gcp_3 = []\n", + "u_gcp_2 = []\n", + "i = 0\n", + "# Ugly way of getting keypoints that are in all 3 images\n", + "for uu in u2[inliers]:\n", + " j = 0\n", + " for vv in u2_p:\n", + " if np.all(uu == vv): # a match between i2 and i2. We can use this\n", + " # to get the real world coordinates computed by correspondences b/t \n", + " # i1 and i2.\n", + " u_gcp_3.append(u3[j]) #camera coordinates to be used for triangulation\n", + " u_gcp_2.append(u2_p[j]) #camera coordinates\n", + " X_gcp.append(real_world_coords[i])\n", + " j+=1\n", + " i+=1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In order to minimize the misfit between the point clouds computed from the two pairs of images, we need control data." + ] + }, + { + "cell_type": "code", + "execution_count": 42, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAADuCAYAAAAOR30qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsnXmYJHWZ5z+RZ1XW1XV19VF9VXVVd9P0fdANgigqHjPiQMuAM7CIDAOCC7ogXjsz4o7aioMXoijKrKsu6zjK44zoeCCH2E1zNWd33feVWWfeR2TsH9W/6Mg7IiuzO6uIz/P4tFRlRkRmRXzjjff3vt9XUhQFExMTE5Ozj+VsH4CJiYmJyTymIJuYmJiUCKYgm5iYmJQIpiCbmJiYlAimIJuYmJiUCKYgm5iYmJQIpiCbmJiYlAimIJuYmJiUCKYgm5iYmJQINoOvN9v6TExMTIwj6XmRGSGbmJiYlAimIJuYmJiUCKYgm5iYmJQIpiCbmJiYlAimIJuYmJiUCKYgm5iYmJQIpiCbmJiYlAimIJuYmJiUCKYgm5iYmJQIpiCbmJiYlAimIJuYmJiUCKYgm5iYmJQIRs2FTEyyoigKsiwDYLVakSRdniomJiaYgmxSIOLxOLIsE4vFCIfD6s8lScJqtar/s1gsWCwWJEkyxdrEJAlTkE0WRDweJxaLEQqFeOmll9i7d68quIoy79YqhFqLJElYLBasVis2m80UahMTTEE2yQNFUVAUhWg0SjweB8BisagCLBDCmk5gFUVRo+rZ2VlGR0dpaWkxhdrkDY0pyCa6ESIai8VUIRZCKURaL1qBlSSJcDiM1WpNEOpIJJLwHm3qQ6Q/TKE2WUqYgmySk2QhFiKoFUJJklSRNkrydnJF1IqiJLxGK9DJeWoTk8WEKcgmGREVE7FYTBXBTIKpzRkbRc97jQj16OgodXV1lJeXZ1xQNDEpRUxBNklBK8QvvvgiO3fuxGLJXrJ+tkQunVDPzc1RU1Oj5rkjkUjCa0QkbQq1SalhCrKJiqIoxGKxhGgzGAwaFqp8Xp9vdJ1tm+luImI/sViMaDSa8DtTqE3ONqYgm6hCLErTMolZsSiGIGfbl/ZfQTahnpycpKmpKaHqw2x6MSkGpiC/gdE2c0DmPO0bgWxC3dfXR0NDQ0IuXVGUrBH1G/V7NFkYpiC/AREVE6LFudACklwFkYszGSEbRXw3yU8MeppebDabKdQmhjAF+Q2CqBMOhUIMDAywbt26khGIUhbkTORKfcTj8YQWcvFai8VCLBajoqLCFGqTFExBXuIk1xDLssz4+DgbNmwo+n7fiCKjpzvx+eefZ8+ePQm/M5teTMAU5CVLpmYOm81W1GhU27ln5D2LLULOB20dt9VqVX+up+lFm/owKz+WLqYgLzFyNXPkI35nItp9IwhyJvLpTjRL9JYmpiAvEdIJcSFK14SAF/NCN0UkPdmEWtv00tHRwcaNGxPax02hXpyYgrzIURQFv99PLBbD6XQWvIbYYrEQj8eLWpf8RklZQGGeBJKFOhAIqOILZtPLYsYU5EWKtpnD7XYTDAZpbW0t+H6EIBs5LpPMFONpQ9RE59P0IoQ5XYmeyZnHFORFRrpmDpvNlrfTWi6MRq99fX0MDw9jsVhwuVxUVlZSUVFBRUWFGsEvdB+LmWKlf7JtM5dQJ6e6fD4fiqJQW1trNr2cYUxBXiRka+YwGsUaQc+2Y7EYAwMD6oW8Z88eFEUhGAzi8/mYnp5maGhI9TxOFmqHw1GUYy9FRMVLKZBJqAOBANFolOrqarPp5QxjCnIJk24yR7oTv5iCnC16jUaj9Pf3Mz4+TnNzM9XV1axdu5Z4PI6iKFRVVVFVVZXwnlgsRiAQwO/3MzU1xcDAAJFIhGAwyIkTJ1SRrqysXJJCLdILpUw8Hk/ISQv0NL2Yk14WhinIJUi2yRzpsFqtauRsZB96LpR0Yh+JROjr68PtdrNmzRoOHjyIxWJhfHw8weshHTabjerqaqqrqxN+/swzz7By5Ur8fj+Tk5P09/cTjUax2WyqSAuhttvthj7rQih0KmUxNMzIspz2Zqin6SXTpJd4PI7T6VTF2hTq9JiCXEJohfj48eNs375d14lrNEIWr9c2J2RCK67hcJi+vj4mJydZu3atKsTJr83nQpMkiZqaGmpqahJ+Ho1G8fv9+P1+3G43fX19RKNR7HZ7glBXVFScUaHOl8UgyPlU1eSqpT5+/DhbtmzBbrer34HZ9JKKKcglQLoaYr/fD+ir0S2mIFssFkKhEENDQ0xPT7N+/Xra2trSXrBi0GkhLyi73c6yZctYtmxZws8jkYgq1OPj42rpn8PhIBQK4Xa7kWWZiooKbLb8T/NCf55i5JALHcXLsqzr3NCDEOp4PI7D4VDPG7PpJT2mIJ9FsjVziMc8PReG0ZSFXgEPBoPMzMwwPT1NW1sbmzdvzrmaf6aqJRwOBw6Hg9raWvVnIt/+6quvIssyo6Oj+P1+ZFnG6XSmRNSFEh0jFFrgixFxF6PuPPlGpLfp5Y0m1KYgnwXSTeZIvgAsFovuSCXfCDkTgUCAnp4efD4f5eXlrF+/nvr6+pzbFZHQ2RA6sX+Hw4HT6WTFihVUVlYC89+3iKh9Ph/Dw8MEAgFVqLUVHy6XK8VnotACWkixK4YgFzJC1qLnOLMJNSTWUnd3d9Pa2rqkhNoU5DOIkckcRqLeQgmyz+ejp6eHYDBIS0sLW7dupaOjQ3fUW6r1xJIk4XQ6cTqd1NXVqT9XFIVwOKymPqampggEAsTjccrKylSBlmW5YFFjMVIgxYhmS60SJN2C4uzsbM7uxK997Wt84hOfWBTrC2AK8hkhHo/j9XqRJAm73a5roc6IIC80ZeH1eunu7iYSidDa2kpdXV1eZkSlKsiZkCSJsrIyysrKEp4AhG+03+/H6/USDod57rnnUBSF8vLyhLSHy+UyJF6LQZCLFSEXg1xNLz/72c/49Kc/fcaPK19MQS4i2maO/v5+amtrWb58ua73ihyyHoxe4EKQ5+bm6O7uRpZlWlpaEqLH5NfqPY7FJMiZkCSJ8vJyysvLqampYWZmhl27dqnNLiKi9ng8BAIBgBShLi8vzzhktdQFudDbLMY5keuc1J6Liyl1YQpygcnUzGG0vVnkkItBNBrl9ddfx26309ramlLBoMWIyArxXkwXgBEkScLlcuFyuWhsbFR/Ho/HE4R6YmKCYDAIgMvlShBqWZYLKnbF+L4LHSEXI+KWZVl39cxiOh9NQS4QuZo5rFZrShtqNvJp9sjF9PQ03d3d+P1+NmzYwNq1a3O+540YIWvRE9FaLBZVcLXE43G1K9Hr9TI2NobP50OWZaLRaIJQl5WV5SUciyFCNiKeeonFYjm3GYvFFk3qRWAK8gLJNJkj+eKy2WyGBNZIyiLX8U1NTdHd3Y3D4WDTpk2MjIxQXl6u6/2ickLva5eaIEP+EZbFYqGyslKt9gCYnJxkenqaFStW4PP5mJ2dZWRkhFAoZMiQSVCMVuxCC3IxhFGPIM/OzqY0GpU6piDnSa7JHMlYrdaUltJsLDRCVhQFj8dDT08P5eXlnHPOOaowiAYOPRh5rVaQF0NHmh6K0TqdTqhhPpIMBAIZDZlE67gwZBI3y1KriEimGCkLPSI/NzeXNR1XipiCbBAhxGNjY1itVmpra3VdEPlUQuTjTwHgdrvp6emhsrKSc889N+VR2kgawshxlFKEXMjjOFOdelarVZch0+DgoCrU4slrampK9flY6PGWek5abDNXhDwzM2NGyEuV5GaOQCCAxWJJW5mQDqOCnE+EPDo6Sn9/P9XV1ezYsSNjWqJYeeFius4ZodCNHIUknxRDJkOmWCzGyMgIU1NTWQ2ZzrbFabEiZDNl8QYkUzOH3W5PsSDMRj6CnFzonun4RIvw9PQ0u3btoqysLOt7jEbIb+Qccim3OttsNsrLy6murqalpUX9eakZMp1NQTZTFkuEdJM5tBdSsSPeXK+Px+OMjIwwMDBAXV0dNTU1bNy4EafTmXPbxRLZpSjIheZMtGIbNWSy2+1qbtrlcpW0WZEgFovhcrmyvsYU5CVAtskcWmw2W1HL2DLlbuPxOMPDwwwMDNDY2MjevXtxOBw8//zzhqJePdG3eK1R8V4Ki3mCUo6QAUN1zdkMmXw+H36/n9HRUYLBIM888wwOhyOh4iNfQ6Zilb3lOpaZmRnWrFlT0P0WG1OQ0T+ZQ4tRgV2ogMuyzNDQEENDQzQ1NbF///6Ex00ji29GRdaI0MdiMXp7ewmFQmolwdlyVisEpW6/udCyN2HIVFdXR11dHZFIhHA4zI4dOxIi6mRDJu3AgGRDpmRisVjONJpR9Ii8WWWxyNDWEHd2dtLY2EhNTY2uC6bYEbKoQ5ZlmcHBQYaHh1mxYgXnnXde2hPRSN2y0RyynkfYWCzG5OQkPp9PdYcLBAIMDw/j9/uJx+Nqe7EQ6UztxUsZRVEKenMqRhOHcErTY8g0ODiYYsikFWojroVG0JNDnpubMxf1FgOidE24eIkoMBqN6o5eiuVBLIjH48zOznLkyBFWrVqVUYjz2X4hXysGnI6OjuJyuWhpaWH16tVEo9EUwx7RXuzz+RLaiysqKgiHw3g8HiorK3M2Q5xJSj1lEY/HC5oOyGWfqseQKdk5Tzixie5Eo4ZM6TAX9ZYA2Zo57HZ7UVMQei9CMTh0dHQUSZI4ePCgbk/kYqUs0kXIQohHRkZobm7mwIEDDA8PZ/yc2Xwg/H4/s7OzTE9PMzw8rNbYikha/LtYLBSzUermQvluT2vI1NDQoP5cURRefvllKisrCQaDeRkypUOvIGvz5YuBN4QgZ5vMIbDZbLoXuqDwXhORSIT+/n4mJiZYs2YN+/bt4+WXX9b9qFfMlIX2tbIsMzAwwPDwMM3NzQk3jHyqLCwWC1VVVarRkdbfVkTTySOatEKdK3+5UEo9h1yMlEUhv08R8DQ0NCTUxRsxZCovL0/5zvR8blOQSwwjhvA2m81QXbGRluJsZBocKtIpRo6nmLXFIpc9NDTE6tWr00buRhYAM+1HYLPZUoaeiskfoiJAm78sLy9XhVrceEuRQntPFFrgC+1GJ7aZfK5kM2QKBoP4fD7VkCkUCgEktI/rOc/8fn/K9kudJSnIomIiV+maFpvNpg4WLfaxSZJEKBSit7c34+BQo63TxYqQYd7A/siRI6xcuZIDBw5kfFQUVRbFyv9qF5rS5aeFUHu9Xl577TVsNluCWU9lZaXqAaGXUs8hF0Pgi9HmrHebuZzzhCFTOBzm2LFjGQ2ZtNtbTCwpQdY2c/z5z3/m/PPPN7RIZyQnnA8WiwW/38/AwACzs7Ns2LAh4+DQfEznC5lDjsfjDA4OMjAwgCRJWYVYUMgI2eh7RX4aIBQK0dzcjMvlSjDrGRwcJBKJYLPZUvLTmT5bqQtysaosCkkhjlFryCQcDPft26caMvn9fmZmZhgaGuKZZ57h+9//PoFAgK9+9ats3bqVXbt2JaxdZGJwcJBrr72WsbExLBYLN954I7fddlvK6/74xz9y++23E41GaWho4PHHH1/Q5xMsCUHO1Mxh5MQ3uqgn0HuBBQIBgsEgL730Eq2trWzZsqWgF2ahUhbxeJyhoSEGBwdZsWIFu3bt4uTJk7pW8hcqqoVOM2Qy69E2Qmg9irUDT0V+utCUuiAXa0htodMq4hjT/Y137NjBoUOHeM973kNtbS2//e1vGR8f55prrsm5bZvNxle+8hV2796N1+tlz549vP3tb+ecc85RXzMzM8OHP/xhfv3rX7N27VomJiYK9tkWtSCLnGKmZg4jJ7/RRT04LWzZTmDt4NCysjJ27NhRlAvdiL1nOkHWdgBqG0+0328uFkvrtN1up7a2NqVjLRwOq0I9OTlJIBBQn5r6+vpUoc7XTB5KpyriTG2vGOhNgdTV1XHdddcZ2vbKlStZuXIlAFVVVWzZsoXh4eEEQf7xj3/M5Zdfrg540DuWTQ+LWpC1Aw6TLxBRBaG3RtNoGZv2PelODu3g0JaWFurr63nppZeK5oaWb8pCK8TLly9P6QA8U14WZ1vMtfW12rKtqakpRkdHKS8vx+v1Mjo6qprJa1MeIj+di1KPkGVZLvnyQj0lbzMzMwuuQe7r6+OFF17gvPPOS/h5R0cH0WiUiy++GK/Xy2233ca11167oH0JFrUgQ+bHb7vdrloR6iEfQU5X+iYGh8ZiMXWCc7bX50LvBZyPx/HQ0BD9/f1phTif7Z6tHHIxEa3FTU1NCT+XZVkty9NaXwpHNe0ik/YcLHVBXgwDTs9EU4jP5+OKK67gq1/9alrb0+eee47f//73BINBDh48yIEDB2hvb897f4JFL8iZMCqwWcvY3G6kxx9H6u5GaWyENWtQdu5M2MfMzAzd3d0AtLS0pK1/NLpwKI5JzwWst8oiHo8zOjqK1+slEAiwb9++rJGdUT/kUhTVhZDp+7darWk9ioX/g8/nU21RZVlW24rFWkJFRUVBhK8YZW+lPuC02F7I0WiUK664gr/5m7/h8ssvT/l9c3MzDQ0N6g33oosu4vjx46YgZ0NEyAvF8vDD2K6/HpIjW6uV5muuYe7v/56OUyvTGzduzHoS5DNXT++qd66UhfBN7uvro76+HpfLpesEyjdlYVQkSjVCNkomRzXRVjw5OakOEgBSRjMZzU8Xo+ytlBtNoLjjmxRF4UMf+hBbtmzhYx/7WNrXXHbZZdx6663EYjEikQhHjx7lox/9qOF9pWPRC3KmkzefFEQKbje2G29ESid0skzTQw+h/PCHNH7rWzj/23/Lubl8LTj15PQypRaShVjYdU5OTuo6BiPisJhzyJkoRIpB21Y8OjpKa2srLpcrpbZWDDu1Wq0p+elM50CpC+hii5D/9Kc/8cMf/pBt27axc+dOAD7/+c8zMDAAwE033cSWLVt45zvfyfbt27FYLNxwww2ce+65xj9IGha9IGci3whZewFK/f2Q5WKUAEmWqbrtNiLvfjfkqHMspql9cspCURTGxsbo7e2lrq6OPXv26DKvXwilKqqlhDaizTTsVNs27na76e3tTTCS1+anSz0nXQwvZFmWcy6gzszMpOT99fCmN71J1zl85513cueddxrefi4WvSAXMkJOrsxQ1q0DPQJjtSL198/nl3Vs38jxGF2oUxSF8fFxenp6qK2tPSNCLFjI1OlSFfOz0RiSqW1c1E/7fD7V1tTv9/PKK68UzNa00BGtnvRCPtvMVTo6NzdXkJzumWbRC3Im7Ha72gOvFyHiqiA3NOD58pepv+02JI1YpFxOsjwv3jkwUissXm9k4nMwGOTIkSMsW7aM3bt3F9wUPBcibWK0KaeUKcaQ03y+m2QjecEzzzxDa2trRltTrVDrsTUt9RQILF2nN1gCglzICNlut6vGNCLvuuxNb8LZ2UnZkSNI3d1IfX1Y//Vf4VTkqtjtyA88kDNdAYUb46RFURQmJibo7u4mHA5z8ODBMy7EgqWYQ4bCdpkVuioil62pGH47NDSky9Z0KeWQF5sXMiwBQc5EvimLkZERPB5PSt5VOXQIIRfy3Xcjvfgic3NzjK9axcaDB3Vvv1BTRhRFwe12093dTXV1NTt27ODll18+a2IMS1OQS90MKBPC1jS5bTwWi6ndiFpbUzGWSUwDKdTYrWIIsp5tmhFyiWFkUU90q42Pj1NXV6dWImSksRHl7W8nNj1NZHRU9zHlU/aWnEMWQtzT00NlZSU7d+6kvLxcHUdlhEKLjWgM6e/vx+12qwtQlZWVb8hxTek420NgbTZbykRqra3p+Pg4Q0NDqq1pclleOm/ibBRrwOlSHN8ES0CQF5Ky0A4OXb58Oc3NzVRXV+tqgdW7Dy0LSVkoioLH46G7u5vKykq2b9+esLCRjztcIY1k4vE44+PjTE1NUV1dTVtbm2qJqZ0SIS5u8T/xXb+RIuRSy69rbU3tdrvq26C1NRXexMFgMMHyUgh1JltT0RRTSPQIstfrTWnaWQwsekHORDazoOTBoaJtuL+/35DAnolBp7FYTBVil8uVIsT5YlSQMwmJtrxu2bJl1NTU0NraSiQSweVyJfgWi7xmcruxw+FQZ+oBBetiKwTFENBCba/YNzBtflproCMsL8XfcWBgIKOtaTGqLMSxZaPQw2TPFItekDP9YdJd0No5cKtXr04ZHCrczfRSzLpiRVEIBAKMjY1RW1vLueeeW9DpB/lMnk520hM57GXLlrFnzx4kSeL48eNZt5MurxmJRHj99deJRqPqFBBFURKisHzM5Zc6xXBm09umn8vWVLSNBwIBdV5i8jTqYlGKT1p6WfSCrAcxOHR8fDzj+CEwPjWkWBHy5OQk3d3dKIrC8uXL2bJli+596CUf0yBxEU1OTtLV1UVFRYWaw4b57zndxeDxSPT3S6xbp9DQkPp7h8NBWVkZTU1Nat5P28UmjMfD4bAahWkflxfTTL1CUuiKDUVRFiRm6WxNX331VVasWIGiKAm2pkDC2K2F2ppqWcyll0tCkDPlH+PxOB0dHbjdbtasWaPOq8tEvuY/Rl6fTQSnpqbo7u7G4XBwzjnnEAwGmZ6e1r19I+RjaD8zM0NnZycOhyNtxJ7u7/DTn9q49dYy7HaIRuG++0IcOpT6HSdfPJm62NI1RyQvPpXyTL1CUuiKjWJUgIi/TfI06uTZeXptTfXchHw+X0r0vlhYEoKcjBgcKkzhcwmxwGgVhFEynUjT09N0dXVht9vZsmWLKkKRSKQk7DplWeall17CYrGwefPmjCd78g3K45G49dYygkGJU70K3HJLGRdf7E8bKesR0Uzm8tqLe25ujtdeey2h1TjXqKZMFGNRr1AshvFNmUrUMs3Oy2ZrWllZqZahZit9m52dXZQLerBEBFlEZsmDQ2tra1m+fLnuk0w0hpwppqen6e7uxmq1phW6fKoyCinIfr+frq4u/H4/5557bk5vgGQ/5P5+CbsdVYwB7Pb5nycLcqFm6i1fvpxwOExzczPl5eXqxa0d1VRWVpZSkrcYH28Xw/gmo3XIuWxNp6amiEQivPjiiwl/SxFNl5eXL8h682yzJAQ5GAzS3d3N7Ows69evVweHut3uolZNCIxGUTMzM3R1dWG1Wtm0aVPGiNOIlwWcLpPTa9eZadvi+/T7/WzcuBFFUfJ6BFy3TiG50CUanf/5mSCTJ0QoFFLTHuPj42opV3JJnt1uL0qtdqFYDG3OhTpGYWsq7Ae2bNmSYGsqjJgefPBBHn/8cWw2G4cPH2bbtm1ceOGFus5fvQNOAY4dO8aBAwd4+OGHOXTo0II/n2BJCPLExAT19fUpg0ONOr7lM1dPRHZ6LrTZ2VkCgQA9PT20t7fnfKwyMpZJvH4hg07D4TA9PT3MzMzQ2tpKY2MjkiQxNDSka7vJ30FDg8J994W45ZbEHHK6dMWZqkPWWmFqW421j8pahzVFUSgrK1MXE0upwaXUp4UICnkT0tYga/+WIj9977338otf/IInnniCtWvX8tRTT2VNs2nRM+AU5s+Vu+66i0svvbRgn0s9hoJv8Sywfv36tMJV0KkhGRD7yNZMMjs7S1dXlzq3bceOHboikWKW1WnTC9FolN7eXjweDxs2bFCfMAQLmQRy6FCMiy/2Z62yEMdzNhfi0j0qK4pCX18f0WiUUCikNriIFEm6BpczSaEX4YoRIRcaPZ1/fr+f9vZ2rr76aq6++mrd29Yz4BTgG9/4BldccQXHjh0z/gFysCQEOROFmhqSjWyCPDc3R1dXF4qiqNNEjh07pvvEzyeHbCRCjkajdHd3MzY2xrp16zhw4EDaC3whs/JgPlLOJMSljCRJ2Gw2ysrK1AsVMje4pFtELGY0LctySadAioFeYyFtQ1I+ZBpwOjw8zM9//nP+8Ic/mIKciaJODclBOtH0er10dXUhyzIbN25M8A0wIrJGBNbI62VZZnZ2luHhYTZs2JCzCsXoceRDoSPkYkbb2RpcRG46XYNLRUUF8Xi8YHnpYkTIS2HAqdfrpaWlJe99ZBtwevvtt3P48OGiPUksCUHOhN1uN9ToIdB7wbjdcPJkNeXlMlVViULc2tqa1m3KyE3CaKogl9gLE6WBgQGcTicbN26kubk553bzEcuzWbVwtrwn0vkVa+ttZ2ZmCIfDHDt2rCANLqVeZVEs681c3hgLsd7MNeD02Wef5aqrrgLA4/Hwq1/9CpvNxvve97689pfMkhDkYk4NycTDD1u46SYbVutGolGJd797gssvH+Jtb1uf1fbPaBrCCHrm6jU2NrJ//35GRkYWvN1CcrZzyMVCW2/b0NDA7Owse/bs0dXgkqt7rdSrLM6m9WaxBpz29vaq//+6667jL/7iLwomxrBEBDkT+Qhy8tSQdLjdcNNNNoJBCZg/OX7+8yZ+/vMmbr5Z5t57MwvumRRkYV4vxjlpbUXz8bIwgjnCKfu29DS4iO41YSqf3OBS6lUWi21aiJ4Bp8VmSQhypgsmn0U9PSJ+4kQIi8VK4tc3fwz332/l/e+XOf/89O8tpiCLbSuKovphVFZWsmvXrpTHPLGop4eFLuotZoolyOnI5K6mHXqqbXABVNvMQjS46J1wbmR7i2laiN4Bp4KHHnrI8D5ysSQEORMLGeOUDr/fT3d3N4FADFk+kHEb73iHg09/WuaGG+SUyU7FjpC9Xi+Dg4M4nU62bduW0aqz2BGyUUo1Qi4k+S7CZWpw6e3tJRqN4vf71Vl6kiSlbXDRQ6FzyHrEsxjb9Hq9ZqdeKZJPhJzOYCgQCJwS4gCtra1s21bPXXfJfPazVlJHnkrEYvDZz1o5fNjKAw/EuPLKeML2jQhystNaJubm5hgeHsZisbBt27achfD5uL0ZxUh0WaqCXMgIuZDubKIkz+VysWLFCvXnmRpcRBSdbYJLoassihUh59pmPB4v+I3gTLE4jzqJbJ7I+TZ6AGpXnc/no7W1lYaGBnVfN9wgc/iwlcyDrSVCIfj7v7fxlrdE1Eg5n7l62QTZ5/PR1dVFLBajqakJh8OhqytpoV192ZAkiUAggMPhKOgj8NngTKYsjJLuvMjU4KItycvU4BKNRkt+US+X8fxCLUTPNktaWAsXAAAgAElEQVRCkKFwEZbNZiMQCPDKK6+oQrx169aUC6mxER54IMaNN9qAOKGQhdRo+bSZTmOjom4/HA7rPp5MVR/BYJCuri4CgQBtbW3U1dUxMjKie9tGI2S9360sy2rBfCwWQ5Zl1fc2W+VAqUbIheRMCHI6tCOaMk1wmZqaYnJykunp6RTzpXwbXM5G5584hxajWRQsIUEuBMFgkPHxcQKBAOecc05aIdZy5ZVx9u2b4ujRcXy+Nu64w3YqYj79nmQznYXM1YN5vwlhpLRx48aEqL1YUa8eT41AIEBHRwfhcJjtK1dS4XYTX7sWpaEhbeVAch1uqUY2iy1CNkJyg0skEmHt2rU4HI6sDS56J7gUY8BprnMkEAgUdLLOmWbJCHKmCEtPDjYUCtHT06O2XArbTj00NVnYtGmOXbviXHZZhAcftPLFL1pxOObF+DvfiSUs7OXrTxGJROjt7WVycpKWlpYUIyWj2zYqyJnSLOIGMTc3R1tbG9X/+Z80vv/9CDeh0H33IR06lFI5kFyHOz09jaIoTE9PU1lZSVVV1Rn1iJA8HqT+fpR161A0RuqFpNATPopV9qanwUXPBBdZllX/4kKg54Y2MzOzaBf0YAkJcibEwl66E0MrxBs2bGDLli1MTU3hdrt1b1+bc25shE98QuZDH5JVM52FVllIkkR/f79qLdre3p41Z16MCDndol4sFqO3txe32336BuHx0PSlLyGFw6oJctktt+C/+OIUkUuuwx0eHkZRFJYtW5Z2CKoQaDGTrZDCZvvpTym79daEm0jslKVioSPkxVo3rG1w0fpiZ2pwKS8vJxKJEI/HqaqqKsh4Jj0Legvp0isFlrwgC8HUCrKwmZyenk6JNgsxJ6+xETVnrOf16ZBlmYGBASYmJli1apWuqSfFTFmIp494PM7AwABDQ0OsXbs2wZBI6u9HsdlAm8e22+d/riPqlCRJFV1ROaBdkPJ6vbjdbgKBQMKIJ/G/fB6PJY+HsltvRQoGc95EFkoxUhaF3F4+VRbZGlxECquzszNrg4te9NYgmxFyCaCnOSQcDtPb28vU1FRam0nIz7LTSAWCHr+JoaEhBgcHWbVqFatXr6a+vl7XhVKslIUkSciyzMjICL29vaxYsYIDBw6kXBzKunVIyd9dNIqybp2ufaQ7nkwLUrIsJ5jMd3d3qwuIoVCImZkZHA4HTqczezNGfz/pxpqIm0gp55CLEXEXYhFOVG84HA7WrFmjjiTL1OCid4KLnpy0Kcgljs1mIxgMMjY2xtTUFOvXr2fTpk0F878weoFlEk1FURgZGaGvr4+mpibOO+88bDabKjR6KEaErCgKXq+X4eFhVqxYwb59+zLndRsb6frkJ2k/fDjh8V9vdJxtwUbN8VZWIvl8SOvWUdPQkNIsEQwGOXnyJH6/n46OjoQFRJH20FYNKOvWkW6sibiJlGodstheoQW50BG3VuD1THDJ1uBSzC69UmHJCHK6EykSiTAzM8PY2BhtbW1ZhVhQ7Ll6yYKvKArj4+P09PRQX1+fInhGxjgVWpDFlGmAxsZGtmzZknO705deyuwHPoBjZAR5zRpSkuh5oOZ4FQVCISgvB0jI9cLpyKy8vJxVq1aptbjaPOfg4KDqAOhyuaiqqmL54cM03HWX4ZuIUUo9hwyFLRfTU/ZmZIJLKBRCkiS6u7tVsXa5XAnfgSnIJUgkEqGvrw+3201VVRUrV65k1apVut5bbGczEQkqioLH46G7u5vq6mp2796d1lbQyBinhaYsTpyAY8csbN3qQ5JOEo/H2bx5M7FYjNHRUV3blSSJeH098ZUrUeLxNJXZmd+n1pB6PFiOH0cBlObm0zlegcFcb7o8ZzweJzQ4SOj11xk95xy6/u//xTo4SHztWsrWrKFyfJzKysq8osZMFRulVvaWTKFrdxdSh5yuwWV0dJRgMEhNTY268BsIBACoqKjgiSeeoKOjg3379hnal55Zej/60Y84fPgwAJWVldx///3s2LEjr8+WjSUjyJIkEY1G6evrY2JignXr1nHw4EHGxsYMNWLku28jRfqxWIxjx45RXl7O9u3bM/pNwPyJqbf9eyER8u23W/n2t8XFs4wPfvBc7r9/vstudnZWd41wvg0e4n22n/6Usptvhkhk/hdWK2RKkRhYMEzG8bOfUZ1UWRG95hoikQher1eNzKanp9W8pEh7ZFuMOlMVG1B4QS50HXgxUipiLSFdg4vD4aCjo4Pnn3+eb33rW6xbt45HHnkk53b1zNLbsGEDjz/+OLW1tTz66KPceOONHD16tGCfTT2Wgm/xLDE7O8sLL7ygCrE4Eex2Oz6fr6j7ttlsulaoZ2dn6ezsJBKJsHv3bnWxIxtWq5VQ5v7slNcamaknePnlKN/+tgNtQ8sPflDBbbdF2LzZmJfFQjrurFNTlN1yC5IQYwBZRtFGx1p0LhimHGOGygr54otxNjTgdDrVoZknTpygqakJi8WCz+dj+s9/xvvcc8y0t6Ns3pyQmy7zerNWbJR6lUUxKPSA00xPkVVVVVx//fW88MIL3HzzzZx33nlMT0/r2q6eWXrna+wbDxw4wNDQ0AI/TXqWjCBXVVWlLQ3Ld4yTkYtH+FNk8m3w+Xx0dnYSj8dpb2/nlVde0SXGYtvFEkPhGHbvvQ5gQ8rvf/ITC5/9bNyQJ8hCImT78PB8RJyMw4EiSSBJKTnkbNGxdPIk1mefRd67F2XTptM/z1FZkZByUBRsNtt8rvlzn8P+wAPqWwLXX8/Ypz6ljsNyHD/Obosl8aKy2dTtFkNACxWBFvpmUQz0LOrNzc2pi4b5eCJnmqWn5cEHH+Rd73qX4W3rYckIcqboMF/HNyNtn5lEPxAI0NXVRSgUoq2tLaVWU88FkCuH7HbDiy/Ob+eUp3ZOFEVRi/gVRaGiIv0Yp+np+e2eiQhZkiRCK1dCus9qseB/8kkkn0+tssjVUVfzmc9QofGrjd54I+F77gGyV1YkpxzqPvUpuP56pJMnsT/wQEJO3PX979N0880sPyX20urVWO+8M2Gz8UiE56emcP3pT9iHh4mvXUt01aqEm/eZ6BLMxVIacJqPEEP2WXqCxx57jAcffJCnnnoqr33kYskIciaKNTVES/LNIBQK0d3djdfrZePGjdTX1yeIr4h6Fzp5+uGHLdxwg03VFrsdDh1az8iIRE0N7NyZ2CmoKAput5uuri7q6+upqKigpaWFq69W+O53U7d/1VXz+zUSIS9kUVSuqyP0rW9RdtNNp3PIdjuhb30LZdMmxBGIfz0eiSeesPDyyxIVFRJ/+ZcxNm1SKOvtpeKhhxLE0/7AA0T+7u/mt9PQQOi++yi75ZaEXC+QknLY8L/+F573vQ/rs8+mPWbbL39JdNMmVVTDn/40zs99bn67skz09ts57wc/wP5f/4VstyPFYpz8+McZffObcTqdND/5JGvvvns+Qo/FUipHzhSFNgIq1oDTYnXq5ZqlB/DSSy9xww038Oijjy54qnUmlowgZ2sMyVeQjb4+EokkdACec845aY9LpDgWIshu97y1ZzSaaGT0k5+s5yc/EccF3//+vB/z9PQ0HR0dVFRUqBUdk5OTAJx/PrztbXF+97vTEdLb3hZXp57kEyEbrkw49b7YoUP4L774dJXFjh0pUaPHI/GNb9j52tccaA/r7rudXHpplG+d92rafViffZbYqWhW7EcbmVqeey4llaHYbFgHB5FbW9Nu03n33diefBLrkSMQj893KdrtEIkQX78e5xe+MP/5ANupm8yWe+5hzXXXEQ6Hqbv7biyhEMLH1XHzzZxsbp6v9DhVf3smHNMWw/gmPU+tsVjMsP+Jnll6AwMDXH755fzwhz+kvb3d0PaNsGQEGdI/LhspGxPkE1UPDQ0RCARyNp6AscW3TDnk/n4pbbpVuzAXi8H119uor3+W2toYW7duzZi7/o//iPH00/C731l529sSR1AVy6ozE0pDA/Ill6T93U9/auOmm8pOPRWkfse/+Y2dS/58Mz3cnfI7ee/elP0klKWlSWVIsRjymjVIHs98flvzdxN7tz72WOKRnNqGpaMjfdlfMIjjH/4B69vfjuRwoDXVlhwOVobDTJ0aTOvz+dThp0Kgq6qqCm66tBgGnOZ6as33vNMzS+/uu+9mcnKSD3/4w8C8Rjyb4alpISwpQU5HPgsVegVZlmX6+/sZGRmhoaEhwdchG0brhdO9dt06BT33jFgM+q7+Ce/Y9xLxd7wD5QMfgMZGTpyAX/963iTG4ZCorFQYGZFwOuf/dbtPpzvO1KJervd5PBK33FKW8FSQZkv0zTXxn7yTv+DXwHyKI3rttQkLe+lQGhoIf/GLOD/+cTXl0PeRj7BsaGheuG229DnuPHD86Efwox+lLGJKsRgVW7fiSqpfDgQC+Hw+dQExHA4TDAbp7OxURTq5ScIIiyVC1rNNo9e8nll63/ve9/je975naLv5sOQFOR9yCXI8HmdwcJChoSFWr15NS0sLkiTpPqGNRsjpXvuHP1hOCbJCukhRS+vcy1h//3usv/893HUX1235M//79fOALXzpSymBHwCbeZ1PvOXPvPMjLdQ02PD1zvGsJb2DnZZ85+/pEeTMTwWp/IK/UgWZigqiH/xgzvfYfvpTnJ/4hJpyiF5yCeu/+lWk++9HisWIXnkl9h/+UHezSya071dkGcXpBKczY5egaCNOdlo7evQo9fX1+Hw++vv7E5oktO54jtnZnIuGiyFCzpUKC4VClJ+qwFmsLClBXognspZMghyPxxkZGaG/v58VK1aofhOig0gvCxVktxtuusmGLEuIJa7V9DPOSmLYOX3Jx7me73Ixf1B/8jqbTonx6RM7+VC+xkf4CN+Ex4DHIIKTN2PhFuf3+L/S1dx11+kBrm43CVaj4m8Qj8cLPjRz3TpFd4D6Pn5++j/i8Zz1ygm1yaewP/ro/Ld0Kvdr/7d/I/zP/4zzM5+Zb+M+hfymN2F96qk00xX1Efrc54jv22eoykK0YafzLfb7/aoznu+732XTl7+MYrNhicWY+OIXka6+OsXAp9ARcjEGnOrxQs5UHbFYWFKCnAkhsHrzbjabLaFUTlEUxsbG6O3tpaGhgf379yeULYnGEL0YbeBIvsn090vYpRg1jDFGM3V4+Bz/k/fwKPdyO1/hf3AT9/MZPk8jngRhOErm+kqATbzOR/hmwnuczHc63he+gf/g7dz92Qae/NzjNDYo/Lv7YmzWOPG4xM17nuYDB5+FizbQcWqlO5uxT67PmUxDg8K3vhXippvK0PaOJNO+fIx3z/0RxVGt25tC6u9P/VnyD+x25PPPx9/djfTEE1gnJoi95S0omzZh/8IXcB4+jFhlPH1LTDcGNxH51DaMkElAtVNAJI+Hiq98Zd6f+lS36vJPfpJn29rwlpUl2GHGYrEz7mNhBD0TZRa7jwUsMUHO5eBmRJCDwaBaJtbd3c2yZcvYs2dPWqN7o4NLjSwapvtM69YpRMMKF/En/o0raKGXC3iSRjx8gc9wB/dQx0zaCO08Trd72oieiqi1v38m47FEsbOePiI4+Hj8i/zVxCNEsRE99VG+c2w/Vxz7BLu+/ixNF74Txw++Rqy2VvUyHhgYUIdrisdq0YqsN/d86FCMiy/2c/y4BVBoblZ49VULf/yjhakpib/92yjr1w8x7jxK9dSU7qhTqaxMbBRJ+wVE1e0pl1+OWOZ03nFHQsMIkgSKwjt4lBfYzU5eYB9H2cpr/JT389f8lKt4GIDvcgNvr99MA9k/e3KTi56INl0DjORwsLWigviePcRiMdV0aWpqikAgwOTkpLqAqJ3aYlSsCy3Iep62Frv1JiwxQc6E0eYQm82G1+vlmWeeoaKigp07d2bNTRXC1N4IVusU//zWX/Ovv93B/+Ea7uQwzQyov69nJuN7t3CSW/k6v+XtfJR7uZnv4CRMiDJA4ij7M77XTpQ+1hPFjgQ4iBDElfB7BzKzVFH/5H/y7MYRIpe/H+cnP8KmTaeL9YWTl9frZWxsDK/XSzQaZWbGTmfnBG1tdtavr8h4A21oULjkktPf36ZNMpdffvq/T548ZXC0IbX7MBOSzwdlZaQbI64AlJWlz++maRhRFIWnbBfx29ilgMTvuJTfcan6+1/xl7zINh7hr5ix1tPeL9HQkFmQkwU/euONhP/5n3MKci5rUZvNxrJly1i2bJlaSbNq1aqEUU2Dg4NEIhHsdnvK1JZs+y+0IOtJgczMzJgR8mLAiGDOzMzQ0dFBJBJh7969ugYm5jMnz2j3IIDX66WjowOr1cp137yAv9l0ATISz7OXITawke6E12eKab7B7QzTxCrGWcY0y5nkOXZyJ1/mJJv5Orfy3/mm+voodmLYuIVv4KeCB7meXbxIhETBFBH0LMuoYY5Wuon++7088u/dfPacj/PNr/qIb2ynv9/OunU1rF59Ot/3r/8a5o47ak/1aUh8/OMnefObR3A6naoQFGoUUDqUdevmI1vtz4C43U7wjjvg7/4ubaSdqWHkv2JvybivKA6+xKcBhXqmEobgJpNO8O0PPADXXJNbkDM0wKT7HCLizjSqSTu1RbispfMsFqm8Qg84fSN4IcMSE+SFNId4vV7V+7e1tZWRkRHd02vziZD1GgbB/MXy0ksvEQ6HaW9vVx/LLP/7X7Bdey1zVDDMZo6wj4v5A01MkM5VQ/vtrGIcgL/mZyjAW/gj/43/w2O8mS5auZ0vcYh/ZxfHqSDILFXcyAP8L/4nq5m34nyQ6/kQ38dOlCh2HuR6GvFgJYaVOFPUMkMdh/g3/uW1j/GjdzxMhHL+peofiUTgS18K88EPRvF4JO68s45w2KJOf7rnns1cd10zVVUh1X1tbGxMNZzXivRCyr0EKeIViRC+805e3LePtvPPzxitJ9c2i+/5Un7D5/jHTHsDoIYp7vnkGA0NazMeV8YOweeeQ9qf+WlGkK4BJh25BDTd4NN0nsXCACgWi1FVVUUwGCzYPD1TkJcIyYt0Wvx+P11dXUQiEdra2li2bBnhcFgtCtdDvpOkcxGJROju7iYYDLJp0yYaGhoSV8avvJLIW95C+S9/ycrvP8HWZ3/HMKupwY9CiBhWXmMzlYSoZ4LlSamMhMdsoBEP7+dnADzBBSxjhiFWIaHwC97H/+CrBHARw4oVmat4mIv5PYOsZz19NOLhddppPxWplxOinFH6WM95HOOz3E0HbfzYexWdbOa225woisKOHXHs9sRsgd0OAwMW9uwpo6ysLMG8PBqNqiItyr20eemxMZnRUSvnnps9FZBMOvGKvvhi1vcomzYRvfFGNaUgvtMLOMI7eJT/ItWEZgXDfIIvcHXZL3Fd9u9Izz2XUSzTCT5AZNcu3Teh5AaYdORTEZPOs1hMABHXVFdXF8FgUF1A1C7uGtmfnhTI3NycOotxsbKkBDnbol6yIAeDQbq7u/H7/arfhPb1xcwJ53p9LBajr6+P8fFxNmzYoBqmpP18jY0o119P+Pzr+cnOOxJSDf+Ha3AR5Jf8JY9xMT1swEUUPzYcxNQ/viieUzT//0Ke5jjn0kUrF/IUb+W3DLCaDQwiY2Gaap5nJ1PUUkmAONDLWtroxorM67QzTR0/5gP8I3dzlP04iNBDC2/h93SyGZC4664yHntshmg0MUcfjZLxUd5ut2eM1n7yE/jMZ9qx2eLEYhY+9akeLr88nLBAlY2U7j0dbeDhe+4hcugQFe98J9pe7t/wHp7iAF9f9UUGarfRGO/j410f403lL0A0SvSaa7BfdFFa72R1/0mCD/M55GhrK5aRkazHZYRC5XzFBBCn00lTU5P6NKddQNQaW2k7ECsrK9Mumov3mxHyEsFut6tF82Li9MzMDK2trTQ2NqZccEYNcgo1V0874LS5uVm1Ex0ZGcl5PJs3w3duvpct99/MZfyCCvzMUsH/5gb+mU/y93ybpzmft/I4LmK4aWCAFezmldTPA0go7ORldvAyElB/6t8YYCVOLXPs5jgTNFCBDz8VVOFljCZ62MAahvgG/52LeJKfcBUn2UI5AVroSajsUCJRPJ4wd9zxOvfcswW7HWIxia99zc+yZVFkWdLVdGO1WolEaviHf6ggHJYIh+fF5QtfaOXSS/uIRCbp7+8nGo2qQzVVH+MCPFIrBw4QvvFGnN/+dsKTx5s4wgWjb8H/yDO8Ksusr/gmAY8HpbKSigsvRAqFck5ACd9zD5G/+7vEKouZmZIeB5Xs1aJdQNTuMxgM4vV6mZ6eVhcQHQ5Hgki7XC5TkBcj2SLkcDhMR0cHHo8n48TpXNspFMllcto658bGRrXhRPt6PRH4vffK/Ne7HHR03Mjb3lbFU09Z8Nxq4+95kAbcfIAfogC1TNJOJ3vwqO8VW7dy+rE7+V/l1O8Fy5ilgllkwEktx9lGA5OEsfEMu/hLfsEn+DK9tFBOgPu5iX7WcIy9SMgoWIkr0HXsRd78ZhtvelOYmZlltLc7WLu2HEk67eMhy7LaDCFJ6UW6v19KsTl2OCSCwSbOOadR/a7FUE2v18vo6GjWvLQRo6TYxz+O89vfTv2FolBxwQXUffrTKDfeSHzdOuxf+lJqRUeWCSjKpk2qMRKUfquznu1pFxC1hMNhNZoWC4iRSITy8nIURVH/TskCbQryIiAWizE2NsbY2BibN2/W7TdRTLSNJB6Ph66uLqqrqzPWORsxSGpvlzn3XC+rVlWxeXMciHHrrTY8NPB1PsrX+SgA1zf/mvcOfZMWOmlmGCcRYjioIpBSEav9b0nz3xIgY6ecKHF8XMITCdGhlzIe4m8ZoplKgrTRya95J/dxG3/DD/GwnNUMIM2u58CBZvx+P3Nzc3i9E7zyihdZlikvL0/JO2pFWvsdrVmjpKvySkh9ZBqqqa0i0OalQ6EQo6OjVFdX53ReUxoaCH/sYzj/5V8SvgfR7dfy+c8zecUVSMEgznvuSa2CiUR0T0ApxnikUhF4p9OpjmoSdHV1YbfbsVqtTExM0NPToy4glpWVcezYsbwEWc88PUVRuO222/jVr36Fy+XioYceYvfu3Xl9tlwsWUGWZZnBwUGGh4dpamqitraW5ub0RuyZMGojqff1VquVcDjMsWPHcDgcuubq5WtGdMMNcS67LMLjj0u89JJEZSW8971x6uvfysaNlxIOz1+EDbjZwYu8ld9wHQ/RyGTakyNZrMuZV0Ar8xG/zOko2kqUBiaZowoPTQRwcQM/4E8c4N94Px5qeQ+/YcWmsYQOM3Vfp0x1vF4vs7Oz6iOtSDkIkXY65xcHa2tjfO1rPj7yERd2u4QsS3zjG0Hq6mQgu9hkqiIQjl7CeU3kPbXRtLZrM3brrTi//nXSOT8pNhu2oSEkp3N+TqAmQlaA8B136G6dLvT0ET0jyIxur9CdetXV1SlDHkKhEBMTE7z44ot0d3fzrne9i2XLlnHVVVdx880359yunnl6jz76KJ2dnXR2dnL06FFuvvnmoszTgyUmyMKzYnh4mIGBAVauXKmOYnn++ecNbSvfqSGZxjgJAoEAHR0dBAIB9u/fr6v33sgYp3TRdGMjHDqkcOhQopx+/vNj3HXXCmIxCQ8NPGG/hAs++Wba7vky2wNP81VuYQfHsTHfAmwHImgFF0LYeJHtvMo2vFTyXn7JBgbwUEuQcmzEWMkomznBDHX0soYLOMIR9vIefkM7Jzn4nvRlX1pTHbF6rigK4XD4VCTtZXx8nGAwiN1ux2azsW6dj0ceacBiaWH9eqivjyPLiSkPITx68tJWq5XVq1err43H4+pNYnJykr6+PjVSEwLd8M1vUnXbbRAOJ0bKsdh8NcUpI6EEysqIXX991uPRUowIudCddcX2xhBPO+vWrePLX/4yR44c4YUXXmBmZoaZmczNUVr0zNN75JFHuPbaa5EkiQMHDjAzM8Po6Kj6vkKypARZURSee+45li1bluA3oSjKGZsakkmQw+Ew3d3dzM3N0dbWRiAQ0G2EYiRlYaTp5JJL3Pz4xx243c2sWrWSvXvnP+uXvgRHuIADnC75asDNo7yVPbzCJNX8kYuoZI7dvMQmXmU3z89HeZTTQStVzGInxkk28G/8NQd5GgcBljPLN7iJn3EVdbi59esbaGzUX6EiSZL6mLp8+XJg/iZ34sQJotEoTU1NBINBAoHnGR62MDd3Ot3hcrmQJEn9LsW/IietZ/HQYrGo0bkgOS89snUr8R/9iDW/+hXrfvQjFLsdSyxGx513sryxEcXh0N2wkQlxcykUhY6QofADTvV6IdfW1hZ0nt7w8DBr1qxR/7u5uZnh4WFTkHMhSRJ79+5NORGK6Ymc6/WxWIze3l7cbjctLS1s2bKlYFUZ6dAj3n6/X43SN29u5L3vXZPw++98J8oNN9gTgjgPjezjZRpw82Ye41PczS/5S77BR1GQOMr5RLBTyzSbeY0b+S61THKMfTTTz9O8iT0cI0IFU6zAYrHwT1+t4UMfyv25TpyAZ5+1snevzObNp38uvtvp6Wk2btyYkG4Qv/f5fMzNzakpByChu0yYHWVaPNSThsqYl774YoZvuYVIRwdT1dWMx+OMHj8+n5o5eJCaI0eo9HiQNmwwPEuv1CPkQqM3OMr3JpBtnl4mB8lisKQEGYw93mdjobXI8XicgYEB9e66kMVEo4Kc6fOLQv25uTna29sJBoNENLZpwjbziisULrwwwssvWxkclHjsMSs/+9n8qeKhkZ9xJT/jSh7mMqao51n2sopBemhlkDUMsobf8h6SM86SBLfe6uev3mrh7l3oiow/9jE73/nO6dP0ppti3HNPhJGREQYGBlizZg0bN25Me4FkKrUSkazH46G3txdZlhPqYYVIj4+PE4vFiEaj6vaFUOv5WzocDhxtbdDWxnIg8PzzbN26lVBovvtwMBTCL0kovb24xscz5qXTEY/Hz7i95dkkV046HA7nPZNxkEcAACAASURBVEUl1zy95uZmBgcH1f8eGhpi1apVee0rF0tOkLNhZJEu3whZOTV6p7e3N8EzeSEYucmkE28x2WR0dDQhSh8dHVW3qyiKGh1KkkRTk4Wmpvl2kQ9+MM5XvhLl8GEb999/Wij+ml+wkmEmWE4bnVzGIxzlAsYr1iPLClddFSUWg8rKGNu3B9i+fQq7fQafz0dvrwW3u5Lq6mpViJIvuBMnOCXGp/9m3/62jQMHnmfHDid79+7NKVzJWCyWtN1lgUCAubk5ZmZm6Ovrw+fzYbPZWL58ObOzs7hcLnXxMB6Pq9+xoihYrdYEwc6EENGampoEV7JseWmtSDudTnU/hY6QCzmUtBgDTnN93nyd3vTM03vve9/LN7/5Ta666iqOHj1KTU1NUdIV8AYSZFFqZnSRTi9Wq5WpqSk6Ojqora1l3759ubvCdN4gjOSQtRGy9uawatUqtdEk+bWyLKvvEbnUZBob4Z57Ytx1V4yXXrKgKNDcHGd4uIHRUT8+X5ytW8/lnzbUMT4eZe3aeNJkEdep/81XusiyjNfrxev1Mjw8jNfrRVEUKitPi/TRo+nzgHNzW2hvL9ypKxYPHQ4Hs7OzWCwW9uzZg8PhUI9xbGyMYDCoNi0I29CysjJVpMXn0n6P2mg60987W15a7H9kZIRwOKy6rgWDQcN10meKYkwLyUW+Nch65um9+93v5le/+hUbN27E5XLxgx/8oKDHrmXJCXKu9uliCLJYdXU6nTmtOgVGqjjmu9CyOLInvTYejzM5OUlHRwfLli1Le3NQFEWddhKLxaiurqampoaysrKs229shEsumRefcDhMPN5NTU2A9vZ2Nepsbs4dzVut1qzphPHxcSoqhoF9Ke+94ILCnraKojA8PMzg4GDKkNry8nJ18RDm0z5CJD0eD4FAQC3Z03aWaSd1C5EWU1Qgd4WHNi+dbv9zc3NMTEwwOjqaYDQvbhRGBbGUzen1MDc3l9e0ED3z9CRJ4r777sv30Ayx5AQ5E/mkIHJVK/h8Pjo7O4nH46xcuRKn06l7ppdRQdYbIQeDQTweD7Iss2PHjpT6Zu14JdGMki4KrKqqUiNVITACWZYZGBhgfHyclpaWtO3n+aBNJyiKQk3NGJddNsgjj5xedPyrvxpmevokx4+71Lrl6urqjB4IuRB2q+KpJtffw+FwUF9fn9C0oF08HBkZURcPtZaUo6OjVFZWqqkhPZ2H2fY/PT2t1k1n8okQi5fiZpEpvaNnGocRimFOn+u7WQpeyLAEBTmbBadRk/pMc/KEm5Xf76etrY26ujr1kVIvRufq5cohh8NhOjs78Xq9uFwu9fFLIIRYXHxCBEQJWXLXWnKdr2gtjsfjTE1NsXr1avbv31+UrsfZ2Vk6OjqoqqrioYda6OkJaqos6lCUA6oHgjBRD4fDas5V3EiyeVSIv2E0GmXr1q26rVbTkWnxcG5ujr6+PmZmZnA4HESjUTo6OhIWD7N1HmYTaa1IZdp/OmvM8vLyBJEWefFSbgpJ9sVIx1Jom4YlKMiZKEQZWzQapbe3F4/HQ2trK1u3blUveJvNht/v1739QpWyCWe4iYkJWlpaaGtr4+WXX054TfKCXa6Lz+Fw0NDQQIOmFMvj8dDR0YHNZqOiooKxsTHcbneCAGaal6cXrW3jli1b1Jzq5s2wefPpzy9JEi6XC5fLpZqoJzeMjIyMpI32y8rKGBwcZGxsTDWXKgZTU1N0dXWxcuVKtm/fruZ7tYuHyZ2H4jt0OBw5Fw9zRY3azkexAKUoijoNZHZ2luHhYcLhsOr1Mj4+npByyZezMS3EFOQSpZARshBk8Yg+MjLCunXr0pawGbXgNHKDSLdtkffs7++nublZPaZYLJZwEWvzlpkW7LIRDAbp7OxUUyDaSDIWi6n5zP7+fnw+X4IQCP+HXCKtTYG0tram+D7rIV3DCCRG+4ODg8zNzak3HOGrvNAbiZZgMMjJkyexWq3s3LkzISev7TzUiqR28U6YHWVbPBRVGfF4XB1OqudGq72Rab8jr9fLiRMnCAaDuN1uAoHAgvyLi2FUpEeQW1paCrbPs8WSE+RM5JtDHhoaor+/n5UrV3LgwIGMJ1oxPZS1r1UURTUkqqurS5mAbaRyIhui6WJqairFL1pgs9lSuqJE9cTc3ByDg4NqPjVZpK1WK4qiMDExoZYIFiMF4nA4cLlcDA4OUlZWxrZt27BaraoAam8k2gqPXEZCyYjSQrfbraax9KBn8U67eGi32/F6vTQ2NqrDYfVUeGTDarVSVlbG+vXr1Z+JvLSogtE21WgXMNPlpfWkGIxgRshLELvdrjvHqygKMzMz6gTeZNFLh9HJ0/nkkOfm5jh58mTGag6RH45EIrzyyitUVVVRU1NDVVWV7gskHo8zMjLC4OAga9asYf/+/YbEPF31hCzL6qKXKHGTZZloNEp5eTmtra3U1dUVXIy1nXxtbW0JN450RkJCpIeGhhIW5rQ3knTCIG6QK1asYN++fQX5HMmLh6FQSG0PX716NaFQiOPHj6vHmKvCI5tIp+vSy5aX9nq9uN1uenp6VEc+8T1VVVUZshzQgx6Bn5ubMwW5FDEyNSQd09PTdHR0qI92mzQetNnQWmrqwYggRyIRZmdnOXnyZMJMPUHygt3BgwdVK0ux6h+Px9UIUESBySf55OQkXV1d1NfX66o4MPJZRTOE6Bb0+/1s2LABWZbVjjlxjNqcbz7HIPyl+/r6aG5uZt++fTlvKtnK8ETrtdfrJR6Pq1Gi0+lkdHQUu92ekp4oFIqiMDg4yMjICBs3bkzI6ycf48TEBD6fL23nod1uV0VapDkANc2ltx4+nSOfNi8tbmZ2u51QKKS+vry8PO+8tN4IOZ/GkFJjyQlyJnINOhVDTiVJUlfd//znP+vefjFSFrFYjJ6eHjweDzabLa1PR6YFu+RutGRx6ejoUAXQ6XQyPT2N0+lk+/btukv3jBCPxxkcHGR0dJQNGzak9fTQRmDj4+N0dXWp4qK9kWR7WvF6vZw8eZKKigq1uSNf0nX1xeNxvF5vQvVEJBLh5MmTCTeSfMvwtIibcF1dHfv27UsbJWbqPBTfY/LioRBIl8uFw+EgFoupN5VoNGrYES9dXrqvrw+73U5ZWRler5eJiYmUuXpGFoH15pDzMRQqNZacIBuNkIPBoDqIMfmx1gj5zNXLFLFrRzmtXbuW8847j6NHjyZ8NqMLduku3FAoREdHh7q6Hg6HOX78eEGiVO1xejweuru7aWpqyigs4hjFxSq8ArTi4na76e7uJhaLqSItjhNQZyRu2rQpIYorJKJ6YsWKFWzbtk2tnggGg8zNzTE9Pc3AwIDhMjwt0WiUrq4uAoFAXiV5kiSp0XG6xUPx5OT3+4lEIqrPsDhmIOEmL7apNy8tyzJVVVUZ67XFAquoStLmpdOdb+LvnQ1TkEsY4dKlJTmCjUQi9PT0MD09nXG2ntF9GsFqtRJKGuEjFrm6u7vTjnISr1lo5YQwPhLR6rZt2xLKqUS6QxulVlRUJESpekTa5/PR0dGBw+HI+5E+k7gI74epqSlOnjxJKBSioqKCxsZGtUqhEFGqIFf1hIgSM/k2Dw8Pq8clvsPkphttqmX9+vVZx4wZRbt4WFNTo9Z5b9iwQV1AFJNSrFZrwsKdSDfozUtnqrIwmpcWxxAOh3N24YVCoaI82Z1plqQgp0MIstZoJ7lNNh3F8gpIXgQUHWMul4vdu3eniJdITSxEiLVVDU1NTezfvz/hwnG7YWDAxtq1VaxefTrCzJRKEJGNEGoh0tFolO7ubrxeb9qc90IR5WORSISZmRmWL1/O+vXr1TK25GYRbSStNejRgzhfJiYmaG9vN1Q9ka4MLxwOq4uHounGarVSXl7O3NwclZWV7N69u6A3E4G2RXzjxo0JNdja3LQoZ0yusNC7eBiNRnUvbGbLS4tJMR6Ph6mpqQSR1ualRfBVap4e+fCGEWSYvxiOHDnC6tWrs5awCYxODTGCWAQUE0RisRhbtmxJedTWGte89tprqueEKHnSy+zsLJ2dnbhcLnbt2pVywf+//2flwx92CL907r8/wvvfP3+RpUsliFpYsZgkUgmSJBEOh1m1ahU7d+407Mamh1AoRGdnJ7FYjHPPPVd9nLXb7RmbRUSFRygUwul0Joh0plSCtnqiUCV5Yl6cEEBZlunq6mJycpK6ujqi0SgvvPBCXvXc2fD5fJw4cYKqqqqcC7aZyhnFk9PExIS6wCkWD0W07/F4VKc8kZIz4ogHqY0/wWCQ9evXqyV/2rx0b28vv/3tb5FlmRdffJGtW7fqXje4/vrr+Y//+A+WL1/OK6+kTl+fnZ3lb//2bxkYGCAWi3HHHXfwwQ9+UNe280Uy2MNeeF+9IhCNRhMcz8bHx+np6SEUCnHhhRfqFonnnnuOrVu36n7Ufvrppzl48KAuoZyamuK1117DarXS1taWsnoujl3k8rSlY16vV62dFdFpdXW1OilZi7b7ra2tLW1u1e2GLVvKCQZPH7fTqfD008EEQ/hsTE5O0tnZqR6LyBUaXZTLRjwep7+/n/Hx8bQVB3rRirTX61U7+sTxORwOBgYGsFgstLe3F6V6AlBz4qtXr6a5uTnhvNH6Y4i/N5CS388VVMiyrNaTb968OS8Dnkxo8/tTU1NMTEwgSVLC31osHorXa9+rd/Hw+PHjbN68Oe1Tw9zcHE8++ST/9E//xIEDB3j11Ve58847ufLKK3Me/xNPPEFlZSXXXnttWkH+/Oc/z+zsLIcPH8btdrNp0ybGxsbyXSjWFT0t6Qh5amqKzs5O9THw+eefNxTt5ls5kW0fIn87ODiIzWbjwIEDaSsnkvPEdrs9JWrRdsr19vbi9/sT8n/i8V2IV6YbxcCABbsdtNYd4TCcf3453/nO6Ug5HSLCt1gs7NixI21ttIishABpRVqISy6RFu8tRLTqdDppbGxM8e+YnZ1lYGCAubk5tUqgv78/oTW8EI/FoVCIkydPYrFY0j6tQHZz/eRKGW3qSHvDEzfJlStXsnfv3qKMaNJGxbt27aKmpiZl8TBd52F5ebkaaED2vHS2srfq6mrOO+88Vq1axUMPPWTo+C+66CL6+vqyfj5hC+vz+airqyvK07KWJSnIohXUarXy/9v78vCmyrT9O933vaQb3ZOWlq22RcRBK7KICzN8jgjXiCDoIMqiH4J1GJERVFBUPkUc/cEA4yCIjgozMjCAoDMCLasINF0pXSld0qRNm/39/YHv60ma5aTJYT33dfWCJifJe9JznvOc57nv+xk8eDDrUlOBBV+RRH/9L2z90WjDpqamBnFxcbjtttugUCjcYk7YurWkfhvV1dXsRK+trUVHRwc7aa0DS3Kyuc/MTUACnQ6YO9cPRUW9sLZ8oJS8zs5Oh+wUblOOu580s2pra2Mj3a1r0r6+vmzclJBcX+BKpkUDfl5eHry8vGAwGFiGypUU99e/g16ML126BJlMZlP96Aj2aHi0dEQvWgaDgc3IS0tLQ0xMjCAmUGq1GgqFgvHW6WfYUx7S75IqD/k0D6k8nJ4btvZDKJXevHnzMGnSJCQkJKCrqwufffaZIN8jFzdlQO7p6UFmZmafZhKtawkVkO1R36jYJDQ0FAUFBfD392cnDeAZ5gTwCyUrIiICo0aNYpkSbXap1WrU1NSwk4Fb7li/3gtPP+2PK2LGXz7b1/dKBh0b+0sJiI5PSk5Ohkwm65fvhD3mBJUK19TUMBtJqVQKqVQqyMnQ29uLiooKSCSSPgHf19e3D3WLfpeUi6zRaFi9lyu7tl5rZ2cnysvLERsb61GJONfcng4koCwNHx8fdHd3o7m5mfGQ3WlwUphMJnYxzsnJsbjY2oMtwypu87ChoYHR4Oj++Pv7o7GxEeHh4Ra2pdx9l0gkgolC9u7di+HDh+Pbb79FdXU1xo0bh9GjR3u07GONmzIgx8fH2wykzsQh1nDXIY5md4QQ5ObmWhy4lGVBa8RcS0xXodFomKiF2+Si8PX17SMVptkfbdIMHNiD//f/wjB7dj4MBglnuysZNHDlwlJZWcnbO9gVUOYEXXtHRwcyMjIQGRnJRhtduHABBoOhT7mjPzU9s9nMXPJcYU/Y+i659d66ujpoNBp20QkKCoJSqYTZbMaQIUOc8mn7C41GA4VCgeDgYIwYMcImXZJbSuA2OLnNQ2eKOjoVJzEx0aZQyRXYax52d3ejvr4eNTU1LHmhCQ39Tn19fVkis3v3bjQ2NvZ7HfawadMmFBcXQyKRIDMzE2lpaVAoFBgxYoTHP4vipgzI9sBXPs3dvj8ZsvUwUVsnu0QigV6vR319PQsurh7cBoMBNTU1UKlULotabGV/er0e77zTgcWLo+DtbYbRKMELL1SgtVWDCxe64e3t7bZ3sCOo1WrmF8ydl+eIg1xbW8uCNLfc4ShIe5o9YaveS3006urqEBQUBJPJhJ9++snlppwzmM1mXLhwAe3t7cjKyrKbKdozMaI0PMo77+npYWOi6Drp+isqKqDX6232CjwFvV6P6upqBAUFYfTo0fDx8bGYOdjR0YGLFy/i5MmT2LFjB/R6Pfz8/LBy5UqPU1STk5Nx4MABjB49Gi0tLSgvLxfcUe6mZFmYTCabgZROhuDrgdvY2AiDwWDhguUI5eXl7ABPT09HXFycw4Zdd3c3Ojs7WRddIpFYBBV79Umq5GtsbERKSgri4+M9eiBe4SN7ISnJiM7OSrS1tSE8PJzR9CgjgQYVd3wKgF+mYff29kIul7ussuMq5WidUq/XIzAw0KIsQ4OKRCIRlD1B5duhoaHIyMhg2Sq3KadWqy08RvqjjFQqlSgvL0dcXBySk5M9VgbhlrhoY1iv1yMiIgJSqdQjNDxrEELYMZ2VleUwuSCE4IsvvsC6detwzz33wM/PD6dPn8bs2bMxefJk3p85bdo0HDp0CG1tbZBKpfjTn/7EErann34aTU1NmDlzJpqbm0EIQXFxMR577LH+7iKvE+SmDMhms9lmJnzhwgX4+/vzHuHd0tKCrq4uZGZmOtyO1lUrKioQFRXFJLXW2zirE3NZE2q1mrEmuMGvp6cHNTU1GDBgAFJSUgSZXWZtzpOYmGixP7RBQ3+4tDH6w0cmbDab0djYiIaGBqSlpUEqlXrswmIdpC9fvgydTofQ0FDExMS4PfrJFmijU6VSITs7m9eFhauMpPVUa9GNNVWQ3sLr9XpkZ2cLlq1SNoi3tzcyMjIslIdcGh5Xedif47GnpwdlZWXsAuboPVpaWvC///u/CA4Oxtq1a/tNfbwGEAOyNerr60EIQXJyMq/3aW9vR2trK7IdkHG5w0QDAgLg4+ODgQN/mQHnbsOOZiutra24dOkSozlFRESw4OduhsoFHZ8UFhaG9PR03pxhLrdXrVY7FWBQSmJ0dDTS0tIEG4rJLU8MHDjQ4hZdrVZbNLvoj6tBmhDCGA4DBw5EYmKiW38PLnOCrpX6OdAmVlpamtufYw9U0dfQ0OCQDULrvXSN1hk//XE0y4+yTrKyshwyJcxmM7744gu8/fbbWLFiBX7961/faMq8WzcgU09ga1y6dAk9PT2860AqlQr19fUYPHhwn+e6u7tZ9iCXyxEUFISmpibo9XqkpqayRh21xexvw06n06GqqgparRYymQxhYWEWDTm1Wt2njMA3Q+WCKyCRy+W8Oud81m4dpH19faHX69n3FhERIciJxWVPOCpP0GYXt9zBlVxzjYHsfY5CoYCvry/kcrlb7nKOoNFocP78eXh5eSE4OBjd3d3MT9qaOeHu5ygUCoSEhCAzM9PlCyVXZk+/Ty6lka5Tr9ejrKwMkZGRSE9Pd1j+uHTpEp5//nmEhYVh7dq1LtMFrxOIAdkabW1trPnBB5S9wB0Yyh1wan1VpyWOjIwMC7es/lKLqIcCn8nOjjJUKre2dbJyP6e/45P47k9tbS1aWloglUpBCEFXVxe0Wq2FK5o7dCzAUs3XH64vYBmkaWCx9sUICQlhcwXlcrlgTmN0fy5fvtzneON6PtjK+F2ht3E/Jzs726M0Mm4jVqVSobW1lRkGRUVF2V2n2WzGjh078O677+K1117DQw89dKNlxVyIAdka1EA7NzeX1/vodDr89NNPKCgo6DNM1Fa9s62tDY2NjZDJZPD19e23ARCt3yYkJGDgwIH9ap5Y+zjQoMLNqCjLIz4+vt+fw2cd1MXLVvPJ3jr7U0bglic82eSi66S0sZaWFrS2tlqIRDxxMbEG5S7TfgGf/bGmt9Hv01H5SKVSQaFQuPQ5/UFXVxfKysoQExODlJQUizFV9OLs7++PQ4cOwdfXF3v37kVSUhLWrl3Lm5Z4HePWDcgAbI5rspXxOoLJZEJpaSmSkpJQV1eHpKQkm4GLlib0ej0uXLgAtVptIQ+mY5Scdc87OztRWVmJ0NBQpKene/z2l56sLS0tqKurAwA2RZobVDxlCKTRaFBeXg4/Pz/IZDLet9POMlTrIM23POEudDodMzXKyspCQECAzYsJ986kP0HaYDCgsrISWq0W2dnZbnOX6UWPG/x6e3vh6+vLxBZyuRzR0dGCZKBcap4zIUlvby/eeustHDhwgCk1pVIp9u7deyNnx8CtHpD1en0fT2S9Xo8zZ86goKDA6etpZnf69GmkpKTYbHA5athxPRzoj9lstsioQkND4eXlxSY7m81myGQywXi+tuhlXDaCSqViNT+uIRDXWpMPuLJqWid2F9a1Xnp7TgiBwWDAwIEDkZSUJEgNl1KyGhoakJGRYcHjtbWttXkRX4c5aoR14cIFpKam2qRNegrt7e0oLy9HZGQk/P390dXVZcFB5hpWubMGmn1LpVKndy1NTU1YuHAhpFIp3n77bVYGuknGM4kB2XrfzGYzSkpKcMcddzh8LXeYqFqtxq9+9SuL5/vbsOPyUGnw0+l0IIQgPj4eCQkJHh1Jz/1c2s1OS0vDgAEDnHpAcy8mlIrlbCYfle5evHgRycnJSEhIECygtLW1obKyEjExMQgJCWHfq16vd0kk4gz0WIiIiEB6enq/2SDWZQSuSo6ukdIy5XK5ILalwJXzgtq9Zmdn97mbsC4jWPt32HMVtIa1vNpRkmE2m7F161Z88MEHWLVqFSZOnHijZ8O2cGsHZK4FJxeHDx/GqFGjbL6GZqo6nY4Zq1tvbz3Drj8HDpd/m5SUxJzZKHXI2meiv1mK9fik5OTkfgcULl+WnqyEEBakvb290dDQgPDwcJfocq7CWXnCmn+sVqst5NZ8yzLUZL+7uxvZ2dkeYZ1YQ6fTMSaPSqWCn59fHzGLq2wZe+Bm3+np6U4vylxYC0Ws/TusRUydnZ1QKBSsB+LocxobG7FgwQIkJiZizZo1N8XkaDsQAzLfgMwdJkq9ielBRD2OAXjEAIgGyOjoaGYAY2vttCNNhRe+vr4unah0fJK/vz8yMzMFmUBhNpuhVCpRU1PD1mh9onpK0eUOe4JrXESDir0gzQ1cQqgguaC38zExMUhLS4OXl5dDr+b+KiO1Wi3Kysrg5+fnseybawxEEwnq1GY2m5nlq70EwGw245NPPsGHH36IN998ExMmTLgZs2IuxIDsLCDTScj19fVISUnpo0gDgJKSEgwfPpw93t9ATAOkr68vZDKZy40n6tfLveXlNrnCw8Ph5+cn+PgkCq58m5txmUwmi1tzW0b6rvoKU19fT7InrIM0zaQNBgMCAgKQlpaGyMhIQTJ9o9HIqJPZ2dlOewZchZy1MtKRKRBXjtxfCiBf0NmGMTExrCZtreaTSCQIDw+HSqXC/PnzkZqaijfffFOQY9TZNJCtW7di9erVbH0ffvghhg0b5vF1cHBrB2TqpGaNI0eOYMSIESxTjY2NRVpamt1hogqFAiqVigU9VwMKNUvp7u72aICkzSNukNZoNDAajYiJiUFiYqJHGRNcUJVdTEwMUlNTnZZBnEnC7ZVlrhZ7gnKk29ramIqTrtV64omrDU4u6EzDmpoat7Nve/J1GqC9vb1RW1uL8PBwp3Jkd2A0GlFZWYne3l4MGjSoj4zbbDazTHrfvn1Yu3Yt2traUFBQgAceeACTJ09GSkqKx9flbBrI4cOHMWjQIERGRuJf//oXli9fjpKSEo+vgwMxINsLyNQaMTMz02YN0rphZzab2YGvUqnQ09MDHx8fiyBtXUIwmUyor6/n3UhzB+3t7aiqqkJ0dDSkUqmFgQ13YjSl3/X35KQ1dkII5HK5Wx4K3LokDdK0LBMaGsosNykdSyhwp2rYozTaanC6OoVbq9VCoVDAx8dHMEWf9dQTOnmbm0m7y5rggjZW+Vxc6urqMH/+fGRkZOD1119HXV0dTp48idtvv523LsBV1NbW4sEHH7QZkLlQKpUYPHiwIBaeHNzaAdna8Y2OGuro6MDgwYNtUpdcadhxMxSVSmVRQjCbzWhtbWVNDaGyE+q3TOfy2QqQzppx4eHhTuu8XDWfkLe+er0ejY2NqKurY8Iabv2Uqg09NUaJelVTTjFfcL0muJRGW65t1K+hubnZJd/l/oAKSbgUMy5rgsrs3aW2UXMjg8GAQYMGOexPmM1mbNq0CRs2bMDbb7+Ne++996rVivkG5DVr1kChUGDDhg1CLkcMyEajkdVUlUolZDIZWlpakJSUZFE68MTEDspbrqyshLe3NzOgt85OPWHqTsc09ZfnS28juXVeav3JLcsAv8yyc0c1yAf2yhO2JOHumAHR2ndTU5Nbg1Jtva/1hY/WpENDQ5GamoqIiAhBLs5Go9GCEeKsJu1OkKYTxvm48128eBHz5s1DdnY2Vq9eLQhTxRH4BOSDBw/imWeewX//+1+hPTJu7YBMA3FTUxNSU1MZJ7a8vBzR0dGIiYnx2Oik3t5eVFVVwWAwWBjz0MYRt87LFYfwyU65oC5c9fX1Huf5Wtd5KUfa398fAwcODy0qdwAAIABJREFURHR0tEdd5ShcZU84k4Q74h6rVCqUl5cjKipKUIc5GiDVajWSk5Mtxj7RTNoRn9sV0LKBuy5zXMMqLv+Yy+qpr68HAGRnZzssuZjNZmzcuBGbNm3Cu+++i6KiomvCoHAWkM+cOYPJkyfjX//6F+RyudDLubUDslKpRHNzc5+mU01NDZua4G4gpv4W7e3tzJjHGayz066uLgsWQnh4uM3shM7Li4yMtNmE9BToFBK1Ws1c8axd5WgW7a53g6fYE/ZUfLQZFxwczPyQ+WSQ7qC1tRVVVVV2A6S1ST0tIVlTBZ0Fab1ej/LycpjNZmRnZwtCa6RBuqmpCW1tbfDx8eljBGXd4L5w4QLmz5+P3NxcrFq1StDv2hkcBeS6ujqMGTMGf/3rX+3qEjyMWzsg2/JEpvW8rq4uDBw4kJfiyBa4gz5tGbi7CqPRaBFMNBoNq50GBASgtbUVXl5ezOZTCPBV2XGzU5VKxTwmuEHaWcOKGp8LyZ6gdyd1dXVoaWlhHOng4GC2Vk+VkIBf9on+nVwto3A5vV1dXQDQJ5P28vKyMJ9yJuN2Fzqdjk1vz8rKgq+vb59mbE9PD9ra2vDll1/C398fJSUl+OCDDzBmzBjBsmJnlDbaF6ipqYHJZMKAAQPw+uuvW0wDefLJJ/H3v/+dMTx8fHxw/PhxQdb7M27tgGzt+EYbdnq9Hk1NTTYzPnsWlVxwM9XU1FTBFGk9PT2oqqpCZ2cn/P39YTabXQ58fNFfU3rAfnZqbVhEZ6O5a43JF9SvOjg4GBkZGUz00R9JuLP9p1xfT9akqfk7t85Pk4zAwEBkZGQgIiJCMIc+enGmQilHOHv2LJYuXQoAiIyMRFVVFaZOnYolS5Z4fG2Ac0rb7t278f7772P37t0oKSnBwoULhaa08YEYkKmfhaM6Mc34aJ1Xp9MhKCjIIvD5+PgwpzgvLy+7jAZPrZueDNzsmxv46FqpCRB3ra4EE2p+T6Xinmq6WNfOaT3aYDAgIiICqampLq+VL6iHglKpRHZ2NsLCHI9sd8ZC4Wan1ujq6oJCoXDb58IZCCGor69HY2MjkpKSAICt1XoOo7vKSKrqo54aju4gTCYTPv74Y2zduhVr167FXXfdxZ4zm82CNYABx+WIOXPmoKioCNOmTQMAZGVl4dChQ2xQ7jUCr4B8006dptJjOpHCXp3Y398fsbGxbPApV8F1+fJlRnqXSCSIj49HXFycIPU64Bf7zbCwMIupy4Dl1GCpVMrWSoPJpUuXGEc4JCSEBWlbJyhVKDY1NSEjI8Op+b2rkEgkCA4ORnBwMKKiopgFZ2JiIrRaLZqbm1FRUQEAFg1Od42VKCMkMTERhYWFvPaJSr1DQ0ORmJgIwLLO39DQwFgoNEiHhISgpaUFarWa9+y8/qK7uxtlZWWIiIjAiBEj+gR9rjKyrq6OKSP5DMvlgu/YJoqqqirMnz8f+fn5+O9//9unlCZkMHaGxsZGizFqSUlJaGxsvNYBmRdu2oCsUCiwaNEiNnAyPz8fhYWFTkeY02ASGBgIvV7P6HIhISEWBz23Ax0eHu4WA0Gr1TKfXWfOWNZrDQkJQUhICBvcyg0m3LVSShsNxrGxsTZPcE+BT3mCG0wuXrzYb5l1b28vG6eVl5fn9gXTy8sL4eHhFtRIutbm5mZGbQwICEBjY2O/JeGOwPUQHjRokN2g7+3tjYiICAvqo73v1V6Q7u3txfnz5xEcHIyCggKnWfGHH36I7du347333uvjhHg9wNZd/43ik3HTliwoDAYDzp07h6NHj+LYsWM4ffo0vLy8kJeXh9tuuw2FhYWQy+UsMJnNZrS1tTmd7GxtANTT0wN/f3+LerSzGi9XcOHJ+qM1jEYjWltbceHCBRiNRhZMrJWGngJlT0ilUpcnUFg3OB2pIrm2okLXpHU6HcrLywFcuQX29/fvtyTcGajAgw5m9US2aWutXl5ekEgk0Gq1yMjIQFxcnMPPqqiowIIFCzBixAisWLFCsLIdH9ysJYubPiBbgxCC7u5unDhxggXpiooKNlbm5MmTWLx4MR544AGXg5R1jZc2t6xrvFxHscTERCQlJQl2i8f1aZDJZEwpZostQbm8dL2uNiyFYk/Y8m3w8vKCVqtFVFQUMjMzBa3p01t5Wt5xBFssBHpBod+tPac+ajrU09PjkUkhjkCHpgYEBDA/aXpBoZl0SEgIW8P69evx+eef4/333xeUJrZnzx4sXLgQJpMJTz75JIqLiy2er6urw4wZM9DS0oKamhp8+eWXuP/++y22+eabb7Bu3TrW1FuwYAFKS0sFWzNPiAGZL7q7uzF79myUl5dj5MiRqKysZFlrfn4+CgoKkJeXh5CQEJeyHW6NlzuNw2AwIDg4GOnp6YiMjBRslDs1suET9K0nh7jig3E12RN6vZ7V9aVSKbuw0AsK9+LnLgOmu7sbCoUCYWFhbhn02LqgWI95UqvVqK6uFtzykxDC/la2hplyLyjl5eVYvHgxent7ER8fj2eeeQZFRUW8p7a7CjpKat++fUhKSkJhYSG2bduGnJwcts3vf/97nDp1Cg0NDWhrawMhBB999JEFpY0Qgnnz5mHPnj0ICgrCpk2beE0JEhhiQOYLQggOHTpkoSgymUwoLy9HSUkJSkpKcOrUKRgMBgwdOpQF6ZycHN4nPWU09Pb2IiEhAUajESqVChqNxqlRkauglK/AwEBkZmb2mx7HZSDQCwpg2Yij+9Wf8oQr4KoUbRms2zKmt2ah8OUdm0wmXLhwAR0dHbyYGv0BvZBQAZPZbEZYWBgiIiL6JQnnA9ogjIyMRHp6usO/ldFoxLp16/Dll1/ijTfegLe3N06cOIGYmBg88cQTHl0XxZEjR7B8+XLs3bsXAPDGG28AAF566SW2zZw5c5Ceno4XX3wRR44cwaJFi3D48GFB1uNhiAHZ0+jp6cGpU6dQWlqK0tJSnD9/HqGhoSxAFxYW9slEuXXO9PR0m4wGqohSqVQWRkXcejSfwM9V2QnlhUwbRu3t7WhqamJBjwYSdxuctkDpZTRT5SvmcGWuIQWtf/OZduEOuPRGWgrpryTcGcxmMytbOWoQUpSVlWH+/Pm4++678corrwhme2qNL774Anv27GEmP5988glKSkqwbt06tk1zczPGjx8PpVIJjUaD/fv3Iz8//6qsz03c2rQ3IRAUFIQ777wTd955J4ArJ1V7eztKS0tRUlKCTz/9lPlMFBQUsEkkL7/8MkaMGGE3I/H19UV0dDS71edyjjs6OlBbW8uMirjZHr2F5ioHU1JSIJfLBQskEokESqUSra2tyMnJQXR0tIUPwuXLl1mD09qlzVVwPSH6Qy+zx0KhggsupS0oKIj5Nzhj4riL3t5elJWVITAwEIWFhewCExAQgICAAKa+4x4HSqUSFy9eZJJw7nHg6GLd1dWFsrIyxMbGoqCgwGlW/H//93/YtWsX1q9fj8LCQs/uuBPwYUds27YNM2fOxKJFi3DkyBFMnz4dZ8+evaY0O09CzJA9DLPZjL1792LJkiWMPqVSqZCTk4OCggIUFBRg6NChLgco2ozklg8kEgmbzhAREcHkrULBFfaEtXk+Fdxwsz17a+XWv901zXEGqrS7ePEiIiMjWZnGU3MNrT+LWnFmZWWxqcquvoe19actf2YvLy8mkBk0aJBT0c/58+cxf/58jBkzBsuWLROMa+8IfEoWubm52LNnD+MZp6en4+jRo4JKyD0EsWRxrXDgwAGEh4ezRoJer8eZM2dYPfqnn36Cn58f8vLyWJDOzMx06Sqv0+lQUVGBnp4eREVFobe3l3X0raXg7gYSrndwf43paY2XG6RtyZYpvczPzw8ymUwQI3cKjUYDhUKB4OBgZGZmWpRCbNEaud7Mrtb6aaYqhNOcdWlGqVSip6cHwcHBiI+Pd9iQNRgMWLt2Lb755husX7/+mja/jEYj5HI5Dhw4wMQ9n376qYWB/cSJE/Hoo49i5syZKCsrw7333ovGxsYbgWcsBuTrFYQQqNVqHDt2DCUlJSgtLUV1dTXi4+NZPbqgoMBmvdlZTZp29Gk92h32gdA8X5qN0rW2t7fDYDAgOjoasbGxjHolhOUnFV1kZWXxrrXznWvIBW0QUim3kKo+k8nESjxZWVkAwNbb3d3NVJzUDzk4OBiLFy/G+PHjsXTpUkGzYmd0NgDYsWMHXnjhBbS0tMDPzw/FxcVYunQpli1bhoKCAkyaNAnnz5/HU089xUpNb775JsaPHy/Yuj0IMSDfSKC3zkePHmVNw46ODsjlchagGxoa0N7ejvvuu8+uYMXW+9qis3Hl1ba8GtwRd7iKjo4OVFRUIC4uDomJiSxIW4st6HrdaRrSz7I3sskVOLP99PLyQlNTExITEwVtEAJXxCQKhcJhM5KqOI8dO4a33noLZWVlSEhIwJ133onf/e53gqnu+NDZKisrMWXKFHz77beIjIzE5cuXb4QyhCsQA/KNDqPRiLKyMvzjH//Axx9/DB8fHwwYMACDBg1iQTorK8tlC0nrzJRmGzTYtbe3w8vLC1lZWYI2t2jZxWQyOfwsLgvFmsdLg7Sz7I7yl/V6PbKzswUVkqjVasaV9vPzAyHEY3MNrWEymVBZWQmNRoOcnByn+/XTTz9hwYIFmDhxIv7whz/AaDTi9OnTCAsLw+DBgz2yJmvwqQ0vWbIEcrkcTz75pCBruA4gsixudPj4+GDIkCHYunUrNm/ejLvvvhsajQYnTpxAaWkp3nzzTZSXlyMyMtKCeudskgjXUIe6h+n1elRVVeHixYsIDAyETqfD+fPnXbIm5Qt6N9DQ0IDMzEyn6jdrFgrwiyqys7MTdXV1fdgH1KWPSy+zxV/2NKjBUWpqKuLi4pjEm9Z4qbGSq3MNbYFm+0lJScjKynI6A3LNmjXYv38/PvroIwwfPhwA4OfnJ7hBuy2zH2s7TGo2deedd8JkMmH58uW47777BF3X9QgxIN8AWLVqFft/SEgI7r77btx9990AfmEk0Ibh5s2b0dzcjLS0NGaolJeXh7CwMLsnLLc88atf/YoFBq41aUNDg11rUlegVquhUCgQGRnplrmRLYoY16WvurqazbQLCgqCXC4XTBUJ/GLm7uXlhfz8fItasjNHufr6ejY5xnquoT2JdUVFBXQ6HYYPH+6UJ/zjjz9i4cKFePDBB/H9998L2ii1BT50NqPRiMrKShw6dAgNDQ0YPXo0zp496/K8yBsdYkC+wSGRSCCVSjFp0iRMmjQJwJWTvbKyEkePHsU333yDFStWQKvVYvDgwSxI5+bmor6+HlVVVYiLi7PJvXVmTVpVVcV7RiCdcUhvrT098JJr+SmVSlFbW8vk73TKRlVVlc1hru4EaS4HnE+2T2HLUY5rAFRTU2Mx146uV6PRoKqqipfEWqfT4a233sLBgwexYcMGDB06tN/76Q6SkpLYPD4AaGhoYLxw7jYjR46Er68v0tLSkJWVhcrKyqvOhb7WEGvItwh0Oh1Onz6No0eP4ujRo/juu+8gkUhw3333YfTo0SgoKEBqaqrLt83cTI/Wo605vCqVChcvXhTcpwG4MkuxoqKCOfVZ74/JZLIwVXLkJucMPT09KCsrs0mb8xRo/VypVKKpqYk1ZCMjIy0k1tbrPX36NBYuXIjf/OY3WLJkiWD8dD7sCaPRiMTERFy+fBmHDx/G008/3YfOtmfPHmzbtg1btmxBW1sb8vLycPr0aaEnQV9NiE09EbYxbdo0DBkyBE888QTOnDnDXO9qa2uRlJTEGob5+fmIiopyOYBSDm9rayuam5tBCLHwaeBjTeoqDAYDKisrodVqXXZK45r/cKXr9uhslA7Y0tKCrKwswW+raQmG1sC51EYqugkICGD89h9//BHHjx/Hn//8ZwwZMkSwdfFhTwBXONh33HEHqqurER0djblz5/ahsxFCsGjRIuzZswfe3t5YunQppk6dKtjarwFu3YDc0dGBRx99FLW1tUhNTcWOHTtsqqK8vb3ZAZucnIxdu3YBuDI5d+rUqejo6MBtt92GTz755KrX3YSEyWRy6NpG69HHjh1DV1eXhcH/0KFDnXbyueY8crkcERERdq1JuZlpf+rJXCtTbiPNHRBC+igN6Xr9/f3R3t6O2NhYZGRkCEoH1Ov1UCgUkEgkyMrKsnsMUvrd5s2bsWPHDrS3tyM4OBhZWVl4++23LRpqngQf9gQAPPfccxg7dizWrFmDNWvWXA/Oa9cCt25AXrJkCaKiolBcXIxVq1ZBqVRi9erVfbajPrDWmDJlCv7nf/4HU6dOxdNPP41hw4Zh7ty5V2Pp1x0MBgPOnj3L+NFnzpxhkzmowb9MJmPBlLInEhISHFp+2rImJYS4NNKJekIEBARAJpMJKhs3Go0oLy9HZ2cnwsLC0Nvb69SoqL/gXmT4TJbWarV44403cOTIEXz00UfIzc1lfYTk5GTBKH58zIBOnTqFlStX4u9//zuKiorEgOwEN2VTb+fOnTh06BAAYMaMGSgqKrIZkG2BEIJvv/0Wn376KXv98uXLb9mA7Ovri7y8POTl5WHu3LkghKCrq4sZ/K9cuZLNAezp6cGQIUOwdOlSJCUlOcxUbRn/0EnLtOZM69Fc6l1AQADzhGhpaWHsCSFB6WWJiYnIyclh+2U9e48yJVwdQcWFTqdDWVkZfH19+8xVtIXjx4/j+eefx6OPPopDhw6xOjblkQsJZ+wJs9mM559/Hps3bxZ0HTcTbsqA3NLSwsa1xMfH4/Llyza302q1bIZYcXExfvOb36C9vR0RERHswKYDEkVcARWQ3HPPPbjnnnsAXMmMVq9ejcmTJ0Ov12Pu3LlobW2FTCZDfn4+8vPzcdtttzkNTjT4cpkHXFFIc3MzNBoNDAYDQkNDkZ6eznv+YH9gMBhQUVEBvV5vk4VijylByxw1NTXQaDTw9fW1oAraahpy+dIymczpOC+tVovXX38dJSUl+Nvf/oZBgwZ5bsd5whl7oqurC2fPnkVRUREA4NKlS5g0aRJ27dp1q2bJTnHDBuSxY8fi0qVLfR5/7bXXeL9HXV0dEhISUFNTgzFjxmDIkCE2zchNJhPGjRvnsCZ9+vRpzJ07F2q1mjUlHn30UQDAzJkz8d1337ETd/PmzYyYfzNg5MiROHbsmEXAMplMUCgUKCkpwddff41ly5bBZDL1Mfh3xkygopCwsDDo9XqYTCbk5OTAYDAwS0qDwdBHCu6uEo6OCEpLS4NUKuWd5fr4+CAqKoqNygIsPTCampr6+F37+/ujqqoKAQEBFnac9lBaWopFixZh2rRpOHjwoCDsDj4oLCxEZWUlG0W2fft2dmcJAOHh4Whra2O/3+IlC164YQPy/v377T4nlUrR3NyM+Ph4NDc3263B0at5eno6ioqKcOrUKTz88MPo7OyE0WiEj48PGhoaoNVqce+997Ka9KpVq/qUQIKCgvDXv/4VMpkMTU1NyM/Px4QJE1gH/q233sJvf/tbD+399QWZTNbnMW9vb+Tm5iI3NxezZs0CcIUmdvLkSZSWlmLt2rUoKytDWFiYhcowMTHRog7LteJMSUmxUKRJpVK2DZWCNzU1sckmXJYE39IBnQvo7e3dR+DRX/j5+fXhc2u1WlaaUalU8PPzg7e3t8UUa+uLSm9vL1auXImTJ09i69atyM7Odntt9uCMzvbOO+9gw4YNMJvNyM3NRWxsLH7/+98jNzfXgj0hwjXclE29xYsXIzo6mgXQjo4OvPnmmxbbKJVKBAUFwd/fH21tbbjjjjuwc+dO5OTk4JFHHsHDDz/MmnpfffUVTp8+zQJ8UVERm0BsD8OGDcMXX3wBmUyGmTNn4sEHH7xpA3J/QQhBW1sbM/gvLS1FQ0MDUlJSUFBQgOTkZHz99dcoLi5Gbm6uS8GRTjahfh3WpQNra1LuMFOh5wICv3CYQ0JCmPWqrSan0WjEkSNHEBsbi48//hjTp0/HwoULPWrfaQ0+dLaDBw/i9ttvR1BQED788EMcOnQIn332mWBruglw67Is2tvbMWXKFNTV1SE5ORmff/45oqKiGDdzw4YNOHz4MObMmQMvLy+YzWY899xzmD17NgCgpqaG0d7y8vKwb98+dHZ2svePjIyEUqm0+/mlpaWYMWMGzp07By8vL8ycORNHjhyBv78/7r33XqxateqaGIDfCDCbzSgvL8fKlSuxf/9+5OTkQKlUWhj8DxkypF/fnzV/V6vVIjAwEIGBgVAqlQgPD4dcLhc02BFCUF9fj6amJmRnZzvkMFOmxJ/+9CecO3cO/v7+CA8Px7Rp0zBv3jzB1siXzkZx6tQpzJs3Dz/88INga7oJcOsG5P7AUU16xowZvAMyzaC3bNmCkSNHsseocurSpUu4/fbbsW/fPovX6XQ6PP744zhx4gSio6Px2WefITU1FcCVE2Ljxo3w9vbGe++9hwkTJnhor69PXLp0CR9//DGWLFmCgIAA6PV6/Pjjj4wfffbsWfj7+1sY/PeHE0z9gy9fvoywsDDodDoL03zqzOYprrFGo0FZWRnCw8ORnp7uNPAfPnwYixcvxowZMzB//nx4e3tDrVajra1NsMnPAD86Gxfz5s1DXFwc/vjHPwq2ppsAty7trT/wRE1arVbjgQcewMqVK1kwBoABAwZg3rx52LdvH6qrq/Hwww/j/PnzFreAGzduRGRkJKqqqrB9+3a8+OKL+Oyzz3D+/Hls374d586dQ1NTE8aOHYuKigpBs7hrjbi4OCxbtoz97ufnh8LCQhQWFmLevHkghEClUjGD/5dffhk1NTVISEhg3OiCggLExMTYrRurVCooFAoMGDAAo0aNYkGXa01qi8oWHh7u8jgnrrIvOzvbqSG+RqPBq6++irNnz2LHjh0WNXq6DiHBxwyI4m9/+xuOHz+O7777TtA13SoQAzIPTJo0CVu2bEFxcTG2bNmCX//613220ev1mDx5Mh5//HE88sgjFs/t3r0bmZmZSEtLw3vvvYe8vDxWr6bYuXMnli9fDgD47W9/ywLPzp07MXXqVPj7+yMtLQ2ZmZkoLS3FHXfcIeg+X8+QSCSIiIjAuHHjMG7cOAC/zKsrKSnBkSNH8N5770GpVFoY/A8fPhxGoxHff/894uLiMHjw4D60OVvWpNT0R6VSoaqqio1z4mNN2t3dzUY3FRYWOsy2CSH44Ycf8OKLL2LWrFlYu3btNbnw8jEDAq4kMa+99hq+++47sQTnIYgBmQeKi4sxZcoUbNy4kdWkAVjUpHfs2IHvv/8e7e3tjAhP6W0vvfQSmpqaMGTIEAwfPhyPPfYYfvzxR4vP4HrG0rl47e3taGxstMi2RV60bUgkEqSkpCAlJQVTpkwBcCWQnjt3DiUlJfjss8/wzDPPoLOzE3fddRfGjh2LoKAgZGVlOQ16Pj4+iIyMtKA62rMm5breNTQ0oK2tDYMGDXI6ukmj0WD58uVQKBT44osvkJGR4f6X0k84o7MBV+rGc+bMwZ49e262yR7XFGJA5oHo6GgcOHCgz+MFBQWszvbYY4/hscces/n6V155BXv37rWoydkSBlhDIpHYfZwvLcnHxwexsbH4y1/+gpSUFAD2PTxuNvj4+GDYsGEYNmwYvL290dbWhtWrV6OpqQmlpaVYvXo1ysvLERUVZUG94+NI58ialAZiyknu7OxkhvTWGTIhBP/5z39QXFyMp556Cu+//76g/hjOjhvay+jt7cWgQYMglUpt0tkWL16M7u5udjd4Mx9HVxNiQL4K4OsHW19fj6SkJBiNRqhUKkRFRdl8rVQqxYwZMyxoSZMmTbIogeTl5eH48eOMlrRkyRJGSwoMDMTp06cF3uvrC9OnT8esWbMgkUggl8uZeoz6RtCG4V/+8hc2QJZr8B8aGupUCh4YGIjm5mbodDqMGDECQUFBTApeV1dnYU164sQJyGQyfPrpp6iursZXX32FtLQ0Qb8Dk8mEZ5991uFxQ3sZzc3N2L59O7766issXboUAPDqq6+y7Rz1XES4AUKIKz8i+gGDwUDS0tJITU0N0el0ZOjQoeTs2bMW26xbt47MmTOHEELItm3byCOPPEIIIeTs2bNk6NChRKvVkpqaGpKWlkb+85//kPHjx7PXvv766+T111+3+/knT54ko0aNYr8HBwd7cvduOphMJlJWVkY2bdpE5s6dS26//XaSl5dHpk+fTtauXUt++OEH0tnZSTQaDftpbGwk3377LTl//jzp7u62eI7709nZSerq6sgTTzxB5HI5SUxMJOPHjycvv/wyMZlMgu7X4cOHnR4348ePJ4cPHyaEXDluo6OjidlsFnRdtwh4xVgxQ74K8PHxwbp16zBhwgSYTCbMmjWrzy3g7NmzMX36dGRmZiIqKgrbt28HAOTm5mLKlClMZvzBBx/g0qVLTmeUcbFx40ZMnDiR/W7Lw0PEL/Dy8kJ2djays7Mxc+ZMAFe+M2rw/8EHH+DcuXMIDg7GkCFD0NDQgLS0NLz88stOJ6FotVqsWLECDQ0N2Lt3L1JTU9HY2IgzZ84IWqoA+M22s9fLcOatIcIzEAPyVcL999+P+++/3+Ix7i1gQEAAaxZaY+nSpey2EYDN7VyhJdny8LiWTaQbAQEBARg5ciRrsBJCsHv3bixcuBC5ubmorq7GmDFjMHDgQAuDfzrHjxCCQ4cO4Q9/+AOeffZZ/PnPf2YBODExkc3aExKEB52NzzYihIMYkG9AuEtLsuXhUVlZ6bDZs3nzZixevJgFjnnz5rGR7Vu2bMHKlSsBAH/84x8xY8YMz+7wdQiJRAI/Pz/s37+fCXjMZjNqa2tx9OhRHDx4EG+99Ra6urogl8tx+fJlBAYG4h//+AeSk5OvyZrd6WWIuErgW9sgYg35ugGfmvTJkydJeno6qaiosHi8o6ODaLVaQgghra2tJDMzk5w5c4akp6eT6upq9n4g5acfAAAGc0lEQVTnzp2zeN2mTZvIs88+22ct7e3tJC0tjbS3t5OOjg6SlpZGOjo6PLzHNy70ej05fvw4eeWVVwSvETuDO70MEW6DV4wVtmglQhBwa9KDBg3ClClTWE2aUo+4tKThw4cz562ysjIUFBRg2LBhuOeee1BcXIzu7m5kZmYiPT0dfn5+mDp1Knbu3MlrLXv37sW4ceMQFRWFyMhIjBs3Dnv27BFs3280+Pr6Ij8/H8uXLxe8RgxcMdMfN24cZDIZxo0bZyHxp8dNUVERwsLC0NzcjGnTpuHhhx9mx83s2bPR3t6OzMxMvPPOO1i1apXgaxbxC0QvCxG8vAs2b96Ml156CbGxsZDL5Xj33XcxcOBArFmzBlqtlvkYrFixAoGBgXjhhReuyb7c6uAzvqyiogISicTCKrasrEzwYa23OHgV4sUMWQSvRs5DDz2E2tpanDlzBmPHjmV1Ykev3bNnD7KyspCZmWkz03r++ecxfPhwDB8+nA1DpfD29mbPib66/LFz5072t5kxYwa+/vrrPtvI5XLmj5GQkIABAwagtbX1qq5ThG2ITT0RvJo9XH/gp556Ci+++CJ7LZ1fSF9bVFTES4Tw7rvvsv+///77OHXqFPv9VhSveAJ8x5dRlJaWQq/Xiyyb6wRihizCwrtAr9dj+/btfbLS5uZm9v9du3axGW4TJkzAv//9byiVSiiVSvz73//GhAkTUFpa6lJdetu2bZg2bZowO3iTYezYsRg8eHCfH751f4rm5mZMnz4dmzZtuir1bRHOIWbIIngJV9577z3s2rWL+TNQA6WoqCi8/PLLKCwsBAAsW7YMUVFRvEQIFBcvXsSFCxcwZswY9pgoXrEPIa1iRVxbiE09EYLg888/72OoVFpaivfff7/PtqtXr0ZDQ4PFc01NTRbilQMHDoi31TzAZ3yZXq/HxIkT8dBDD+G55567Riu95SA29URcO/AVrwDA9u3b+5QrbIlXZs2ahQEDBmDw4ME234cQggULFiAzMxNDhw7FyZMn2XNbtmyBTCaDTCbDli1b3N29qwpHVDZrPPvss3jllVcQERGBffv2MYHP8ePHmZCHWsVSe9jhw4eL9frrBXwJy0QUhohwAXxECIQQolAoSEpKioWBjS3xyrlz58h3331HTpw4QXJzc21+5jfffEPuu+8+YjabyZEjR8iIESMIITe+eGXx4sXkjTfeIIQQ8sYbb5AlS5bY3XbBggVk2rRpNkU8Iq4pRGGIiGsHPuIV4Eozb+rUqRY0O1vilZycHNx1110OZbw7d+7E448/DolEgpEjR6KzsxPNzc03vHiFD5UNAE6cOIGWlhaMHz/+ai5PhAchNvVECAZnhkoA2NgqLkaNGoWffvrJ5c+z1UhsbGy0+/iNAj5UNrPZjEWLFuGTTz6xOUxBxI0BMSCLuGlAXJy6AgCzZs3CP//5TwwYMABnz57ts93WrVuZ0i0kJAQffvghhg0bBgBITU1FaGgovL294ePjg+PHj/d77Y6mnvPB+vXrcf/991tceETceHCVZSFCxDWFRCJJBfBPQkifzp5EIvkIwCFCyLaffy8HUER/CCFzrLeTSCR3AegG8Fc77zkKQBkhRCmRSCYCWE4Iuf3n52oBFBBC2jy9n1ZrKP95/c0SiST+57VnWW2zFcBoAGYAIQD8AKwnhBT3eUMR1y3EGrKImwm7ADwuuYKRAFSEkGYAewGMl0gkkRKJJBLA+J8fAyHkewAd9t6QEHKYEEJpDUcBJAm6B7axCwD1NJ0BoI8ChBDyO0JIMiEkFcALuHKBEYPxDQaxZCHihoFEItmGK9lujEQiaQDwCgBfACCE/BnAbgD3A6gC0APgiZ+f65BIJCsAHPv5rV4lhNgNwg4wG8C/OL8TAP+WSCQEwEeEkI/78Z58sArADolEMhtAHYBHAEAikRQAeJoQ8qRAnyviKkMsWYi45eGoDMLZ5h4A6wH8ihDS/vNjCYSQJolEMgDAPgDzf864RYjoF8SShQgRTiCRSIYC2ADg1zQYAwAhpOnnfy8D+ArAiGuzQhE3C8SALEKEA0gkkmQAXwKYTgip4DweLJFIQun/caUu3ZemIUKECxBryCJuafCoSy8DEA1g/c9UOSMhpACAFMBXPz/mA+BTQsiNozYRcV1CrCGLECFCxHUCsWQhQoQIEdcJxIAsQoQIEdcJxIAsQoQIEdcJ/j9Qh+jYPu9NTQAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "X_gcp = np.vstack(X_gcp)\n", + "u_gcp_2 = np.vstack(u_gcp_2)\n", + "u_gcp_3 = np.vstack(u_gcp_3)\n", + "X_gcp_h = np.ones((X_gcp.shape[0], X_gcp.shape[1]+1))\n", + "X_gcp_h[:, :X_gcp.shape[1]] = X_gcp\n", + "initial_trans = initial_pose_guess[:, 3]\n", + "initial_rot = cv2.Rodrigues(initial_pose_guess[:, :3])[0]\n", + "c_r = Camera_Rodrigues(initial_rot, initial_trans, P_2c) # There is still\n", + "projection = c_r.minimize(X_gcp, u_gcp_2, u_gcp_3)\n", + "# a slight tilt in the final point cloud\n", + "real_world_coords_2 = []\n", + "for xx1, xx2 in zip(u2_p[inliers2, :2], u3[inliers2, :2]):\n", + " a = triangulate(P_2c, projection, xx1, xx2) # x, y, z, w\n", + " real_world_coords_2.append(a[:3] / a[3])\n", + "from mpl_toolkits.mplot3d import Axes3D # noqa: F401 unused import\n", + "fig = plt.figure()\n", + "ax = fig.add_subplot(111, projection='3d')\n", + "for aa in real_world_coords:\n", + " ax.scatter(aa[0], aa[1], aa[2], c='red')\n", + "for aa in real_world_coords_2:\n", + " ax.scatter(aa[0], aa[1], aa[2], c='blue')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There's still a minor rotation offset in the final point cloud, though they are very close. Maybe I didn't understand the problem fully. All I tried to do was minimize the difference between the point clouds computed from the two sets of images. To do this, I fit one of the projection matrices in the triangulate function. " + ] } ], "metadata": { @@ -55,7 +229,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.5" + "version": "3.7.2" } }, "nbformat": 4, diff --git a/pictures/DSC03919.JPG b/pictures/DSC03919.JPG new file mode 100644 index 0000000..bff1860 Binary files /dev/null and b/pictures/DSC03919.JPG differ diff --git a/pictures/DSC03920.JPG b/pictures/DSC03920.JPG new file mode 100644 index 0000000..ceeb844 Binary files /dev/null and b/pictures/DSC03920.JPG differ diff --git a/pictures/DSC03921.JPG b/pictures/DSC03921.JPG new file mode 100644 index 0000000..67a3bb8 Binary files /dev/null and b/pictures/DSC03921.JPG differ diff --git a/structure_from_motion.py b/structure_from_motion.py new file mode 100644 index 0000000..0a81e18 --- /dev/null +++ b/structure_from_motion.py @@ -0,0 +1,253 @@ +import numpy as np +import matplotlib.pyplot as plt +import cv2 +import piexif +from scipy.optimize import leastsq +from glob import glob + + +def compute_sift_and_match(i1, i2): + h, w, d = i1.shape + sift = cv2.xfeatures2d.SIFT_create() + + kp1,des1 = sift.detectAndCompute(i1, None) + kp2,des2 = sift.detectAndCompute(i2, None) + + bf = cv2.BFMatcher() + matches = bf.knnMatch(des1, des2, k=2) + + # Apply ratio test + good = [] + for i, (m, n) in enumerate(matches): + if m.distance < 0.7*n.distance: + good.append(m) + + u1 = [] + u2 = [] + + for m in good: + u1.append(kp1[m.queryIdx].pt) + u2.append(kp2[m.trainIdx].pt) + + u1 = np.array(u1) + u2 = np.array(u2) + u1 = np.c_[u1,np.ones(u1.shape[0])] + u2 = np.c_[u2,np.ones(u2.shape[0])] + return u1, u2 + +def u_to_x(u, im, f_im): + ''' + Converts camera coordinates to + generalized image coordinates. + ''' + h, w, d = im.shape + exif = piexif.load(f_im) + f = exif['Exif'][piexif.ExifIFD.FocalLengthIn35mmFilm]/36*w + cu = w // 2 + cv = h // 2 + K_cam = np.array([[f,0,cu],[0,f,cv],[0,0,1]]) + K_inv = np.linalg.inv(K_cam) + x = u @ K_inv.T + return x, K_cam + +def triangulate(P0,P1,x1,x2): + # P0,P1: projection matrices for each of two cameras/images + # x1,x1: corresponding points in each of two images (If using P that has been scaled by K, then use camera + # coordinates, otherwise use generalized coordinates) + A = np.array([[P0[2,0]*x1[0] - P0[0,0], P0[2,1]*x1[0] - P0[0,1], P0[2,2]*x1[0] - P0[0,2], P0[2,3]*x1[0] - P0[0,3]], + [P0[2,0]*x1[1] - P0[1,0], P0[2,1]*x1[1] - P0[1,1], P0[2,2]*x1[1] - P0[1,2], P0[2,3]*x1[1] - P0[1,3]], + [P1[2,0]*x2[0] - P1[0,0], P1[2,1]*x2[0] - P1[0,1], P1[2,2]*x2[0] - P1[0,2], P1[2,3]*x2[0] - P1[0,3]], + [P1[2,0]*x2[1] - P1[1,0], P1[2,1]*x2[1] - P1[1,1], P1[2,2]*x2[1] - P1[1,2], P1[2,3]*x2[1] - P1[1,3]]]) + u,s,vt = np.linalg.svd(A) + return vt[-1] + + +class Camera_2(object): + + ''' Fit the entire projection matrix.''' + + def __init__(self, cam_mat): + self.cam_mat = cam_mat + + def predict(self, uv): + return self.cam_mat @ uv.T + + def fit_func(self, x, cam_mat): + cam_mat = np.reshape(cam_mat, (3, 4)) + self.cam_mat = cam_mat + return self.predict(x) + + def err_func(self, cam_mat, uv, real_world_coords, fit_func, err): + out = fit_func(real_world_coords, cam_mat) + return (uv - out.T).ravel() + + def estimate_pose(self, X_gcp, u_gcp): + err = np.ones(X_gcp.shape) + out = leastsq(self.err_func, self.cam_mat.flatten(), args=(u_gcp, X_gcp, self.fit_func, + err), full_output=1) + self.cam_mat = np.reshape(out[0], (3, 4)) + return out[0] + + +class Camera(object): + ''' Fit just the translation vector. ''' + + def __init__(self, cam_mat): + self.rot_mat = cam_mat[:, :3] + self.trans_vect = cam_mat[:, 3] + self.trans_vect = np.expand_dims(self.trans_vect, 1) + + def predict(self, uv): + return np.hstack((self.rot_mat, self.trans_vect)) @ uv.T + + def fit_func(self, x, trans_vect): + self.trans_vect = np.expand_dims(trans_vect, 1) + return self.predict(x) + + def err_func(self, trans_vect, uv, real_world_coords, fit_func, err): + out = fit_func(real_world_coords, trans_vect) + return (uv - out.T).ravel() + + def estimate_pose(self, X_gcp, u_gcp): + err = np.ones(X_gcp.shape) + out = leastsq(self.err_func, self.trans_vect, args=(u_gcp, X_gcp, self.fit_func, err), full_output=1) + self.trans_vect = out[0] + self.trans_vect = np.expand_dims(self.trans_vect, 1) + self.cam_mat = np.hstack((self.rot_mat, self.trans_vect)) + return out[0] + + +class Camera_Rodrigues(object): + + ''' Fit the projection matrix, constraining + it to 6 free parameters. ''' + + def __init__(self, rot_vect, trans_vect, P_2c): + self.trans_vect = np.expand_dims(trans_vect, 1) + self.rot_vect = rot_vect + self.params = np.hstack((self.trans_vect, self.rot_vect)) + self.P_2c = P_2c + + def fit_func(self, x, p): + u2 = x[0] + u3 = x[1] + trans = np.expand_dims(p[3:], 1) + rot = p[:3] + rot_mat = cv2.Rodrigues(rot)[0] + cam_mat = np.hstack((rot_mat, trans)) + real_world_coords = [] + for xx1, xx2 in zip(u2[:, :2], u3[:, :2]): + a = triangulate(self.P_2c, cam_mat, xx1, xx2) # x, y, z, w + real_world_coords.append(a[:3] / a[3]) + real_world_coords = np.asarray(real_world_coords) + return real_world_coords + + def err_func(self, p, y, u2, u3, fit_func, err): + ''' Gets the error between point clouds. ''' + # use p as a projection matrix to triangulate points. + x = [u2, u3] + out = y - fit_func(x, p) # this line returns the error between the + # two point clouds - one estimated from i1 and i2 and the + # other estimated from i2 and i3. + return out.ravel() + + def minimize(self, X_gcp, u2, u3): + ''' + X_gcp are the real-world coordinates corresponding to image + coordinates u2, u3 + ''' + err = np.ones(X_gcp.shape) + out = leastsq(self.err_func, self.params, args=(X_gcp, u2, u3, self.fit_func, err), full_output=1) + p = out[0] + trans = np.expand_dims(p[3:], 1) + rot_mat = p[:3] + rot_mat = cv2.Rodrigues(rot_mat)[0] + proj_mat = np.hstack((rot_mat, trans)) + return proj_mat + + +if __name__ == '__main__': + files = ['pictures/falcon/DSC03919.JPG', + 'pictures/falcon/DSC03920.JPG','pictures/falcon/DSC03921.JPG'] + i = 0 + f1 = files[i] + f2 = files[i+1] + f3 = files[i+2] + i1 = cv2.imread(f1) + i2 = cv2.imread(f2) + i3 = cv2.imread(f3) + # First two images + u1, u2 = compute_sift_and_match(i1, i2) + x1, k_cam1 = u_to_x(u1, i1, f1) + x2, k_cam2 = u_to_x(u2, i2, f2) + E, inliers = cv2.findEssentialMat(x1[:,:2],x2[:,:2],np.eye(3),method=cv2.RANSAC,threshold=1e-3) + inliers = inliers.ravel().astype(bool) + n_in, R, t, _ = cv2.recoverPose(E, x1[inliers,:2], x2[inliers,:2]) + P_1 = np.array([[1,0,0,0], + [0,1,0,0], + [0,0,1,0]]) + P_2 = np.hstack((R,t)) + P_1c = k_cam1 @ P_1 + P_2c = k_cam2 @ P_2 + real_world_coords = [] + for xx1, xx2 in zip(u1[inliers, :2], u2[inliers, :2]): + a = triangulate(P_1c, P_2c, xx1, xx2) # x, y, z, w + real_world_coords.append(a[:3] / a[3]) + + u2_p, u3 = compute_sift_and_match(i2, i3) # u2 prime and u3 + x2_p, k_cam2_p = u_to_x(u2_p, i2, f2) + x3, k_cam3 = u_to_x(u3, i3, f3) + E, inliers2 = cv2.findEssentialMat(x2_p[:,:2],x3[:,:2],np.eye(3), method=cv2.RANSAC,threshold=1e-3) + inliers2 = inliers2.ravel().astype(bool) + n_in, R, t, _ = cv2.recoverPose(E, x2_p[inliers2,:2], x3[inliers2,:2]) + P_2_p = np.array([[1,0,0,0], + [0,1,0,0], + [0,0,1,0]]) + P_3 = np.hstack((R, t)) + P_2_pc = k_cam2_p @ P_2_p + P_3c = k_cam3 @ P_3 + append_row = np.ones((4)) + P_3c = np.vstack((P_3c, append_row)) + initial_pose_guess = P_2c @ P_3c.T + X_gcp = [] + u_gcp_3 = [] + u_gcp_2 = [] + i = 0 + # Ugly way of getting keypoints that are in all 3 images + for uu in u2[inliers]: + j = 0 + for vv in u2_p: + if np.all(uu == vv): # a match between i2 and i2. We can use this + # to get the real world coordinates computed by correspondences b/t + # i1 and i2. + u_gcp_3.append(u3[j]) #camera coordinates to be used for triangulation + u_gcp_2.append(u2_p[j]) #camera coordinates + X_gcp.append(real_world_coords[i]) + j+=1 + i+=1 + + X_gcp = np.vstack(X_gcp) + u_gcp_2 = np.vstack(u_gcp_2) + u_gcp_3 = np.vstack(u_gcp_3) + c1 = Camera_2(initial_pose_guess) + X_gcp_h = np.ones((X_gcp.shape[0], X_gcp.shape[1]+1)) + X_gcp_h[:, :X_gcp.shape[1]] = X_gcp + c1.estimate_pose(X_gcp_h, u_gcp_3) + initial_pose_guess = c1.cam_mat + initial_trans = initial_pose_guess[:, 3] + initial_rot = cv2.Rodrigues(initial_pose_guess[:, :3])[0] + c_r = Camera_Rodrigues(initial_rot, initial_trans, P_2c) # There is still + projection = c_r.minimize(X_gcp, u_gcp_2, u_gcp_3) + # a slight tilt in the final point cloud + real_world_coords_2 = [] + for xx1, xx2 in zip(u2_p[inliers2, :2], u3[inliers2, :2]): + a = triangulate(P_2c, projection, xx1, xx2) # x, y, z, w + real_world_coords_2.append(a[:3] / a[3]) + from mpl_toolkits.mplot3d import Axes3D # noqa: F401 unused import + fig = plt.figure() + ax = fig.add_subplot(111, projection='3d') + for aa in real_world_coords: + ax.scatter(aa[0], aa[1], aa[2], c='red') + for aa in real_world_coords_2: + ax.scatter(aa[0], aa[1], aa[2], c='green') + plt.show()